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Minimum Weight Analysis Based on 
Structural Reliability’ 


HARRY H. HILTON* ann MORRIS FEIGEN** 
Unwersity of Illinots and Hughes Aircraft Company 


Summary 


An analytical investigation is presented for the proportioning 
of probabilities of failure among structural components in terms 
of a preassigned probability of failure of the entire structure, 
such that the total structural weight is a minimum. 

The solution is expressed in terms of general loads of known 
probability of occurrence and normal distributions of failure 
stresses, which are assumed to be known from prior tests for 
given types of structures, loads, and materials. The types of 
failures considered include simple tension or compression, ten- 
sion or compression creep, elastic and creep buckling of cylinders 
and plates, bending and torsion of round tubes, and their com- 
binations. It is shown that general relations can be derived for 
these loadings and, consequently, the analysis is in terms of 
generalized loads applicable to the above cases. 

A solution is obtained for the structure consisting of many 
elements, where each element may have different material prop- 
erties, standard deviation and loading. From the exact solution 
empirical relations between the ratio of probabilities of failure 
and load ratio for pairs of elements are deduced. Similar empiri- 
cal relations were also found to exist between weight ratios and 
probability of failure ratios 

It is shown that for a given total structural probability of failure 
of a multicomponent structure, as any component becomes heav- 
ier it should be assigned a higher probability of failure in order 
to achieve minimum weight for the entire structure. 

The relationship of the margin of safety to probability of failure 
is discussed. Results of numerical computations determining 
the aforementioned empirical relations, the margin of safety, and 


the weight saved are presented. 


Received July 20, 1959. 

} The analysis described in this paper was performed under the 
sponsorship of the Systems Development Laboratories, Hughes 
Aircraft Company. The authors are grateful to HAC for per- 
mission to publish this article. Portions of this paper were pre- 
sented at the Third Symposium on High Speed Aerodynamics 
and Structures, San Diego, March 27, 1958. 

* Professor of Aeronautical Engineering, University of Illinois, 
and Consultant, Hughes Aircraft Company. 

** Senior Staff Engineer, Hughes Aircraft Company. Now 
Senior Staff Engineer, Space Technology Laboratories, Inc. 


Symbols 


cross-sectional area 

[Crm? + Cre? — Cim2Crm*Gm?]"2 

outer radius of round tube 

ratio of standard deviation of applied stress to mean 
applied stress 

ratio of standard deviation of failure stress to mean 
failure stress 

mean applied stress 

applied stress 

design stress 

mean failure stress 

failure stress 

minimum guaranteed value of failure stress 

difference between mean failure and mean applied 
stresses 

difference between failure and applied stresses 

nondimensional form of gm 

nondimensional form of g,»,’ 

constant 

length 

constant exponent in Eq. (17) 

probability of failure of mth structural element 

probability of failure of the entire structure 

see Eq. (46) 

number of structural elements 

thickness 

weight of entire structure 

weight of the mth structural element 

weight saved 

orm" [orm 

OFm f oFm 

percent of material 

see Eq. (21) 

mean applied load 

applied load 

Lagrangian multipler 

density 

standard deviation of applied stress 

standard deviation of failure stress 

standard deviation of difference between failure and 
applied stress 

difference between design stress and mean applied 


stress 
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or,,” = difference between minimum guaranteed value of 
failure stress and mean failure stress 
db». = nondimensional statistical distribution 
€n = statistical distribution 
Q = weight parameter, see Eq. (48) 


Barred quantities refer to conventional analysis based on 
margin of safety 
Subscripts m identify the mth structural element. 


Introduction 


i be ANALYTICAL AND EXPERIMENTAL EFFORTS in 
aircraft structural research and design have always 
been directed toward the development of a_ high- 
strength, low-weight airframe. However, not until the 
last fifteen years has the ever-present problem of struc- 
tural weight reduction received increasing rational and 
systematic attention. 

Two main approaches have been used to achieve this 
goal. One has been the initiation of optimum struc- 
tural analysis practices, which have sought to produce 
structures of maximum strength-to-weight ratios by 
proper selection of materials and configurations." ” 

The second approach has dealt with the re-examina- 
tion of failure and safety factor criteria and their 
interpretation in terms of structural reliability. This 
necessitates the continuing collection of statistical data 
for applied loads, failure stresses and environmental 
conditions. Freudenthal* has discussed in detail the 
probability of failure criteria, and English’ ° has ex- 
tended these to include statistical considerations of 
tolerances and cost. 

The structural reliability concept is particularly well 
suited to missiles because of their limited-purpose 
missions, relatively large number of missiles per mission, 
short life expectancy, and lack of human life involve- 
ment except during captive flights and surface launches. 
In other words, instead of overdesigning the missile 
to meet specifically the one worst condition which 
infrequently, considerable weight 


generally occurs 


can be saved by allowing probable structural failures 
in one out of X flights. 

Generally, for every pound of structure eliminated, 
the overall weight of an air-to-air missile can be re- 
duced by approximately 2 to 3 pounds. 


For an air- 
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plane, the equivalent ratio is 5 or 10 to 1. Conse- 
quently, an airplane equipped with missiles can have 
a saving of 10 to 50 pounds of total weight for every 
pound of structural missile weight eliminated. 

One question that still needs to be answered is how 
much one is willing to spend to eliminate each missile 
structural pound. Obviously, the optimum design is 
not necessarily the cheapest one, since in order to 
achieve it one may have to use more costly materials 
and/or fabrication Finally, the growth 
factor of the missile may also enter into the design con- 
siderations. The optimum structure may not lend 
itself to future expected missions and one may choose 


processes. 


to overdesign for the present in order to satisfy pro- 
jected requirements. 

In the present paper, the concepts of minimum 
weight analysis and structural reliability are combined 
to answer the question of how to assign probabilities 
of failure to the various structural components such 
that the total structural weight is a minimum and 
the probability of failure of the entire structure is equal 
to a preassigned value. The results presented herein 
should help the analyst and designer proportion struc- 
tural members and obtain a minimum weight structure 
whose reliability meets predetermined standards. 


Analysis 


Consider a structure consisting of r elements, each 
of which may have different loading and environmental 
conditions, material properties, and geometry. It will 
be assumed that the mode of failure and the statistical 
distribution of failure stresses and of applied loads (Fig. 
1) are known for each element and that the latter have 
been determined prior to the present analysis. Failure 
stresses rather than failure strengths have been chosen 
as representative quantities, because then the statistical 
distribution (¢,,), standard deviations (o,,,), and mean 
failure stresses (/’,,) are fixed material properties and 
known constants for any element where the allowable 
stress is independent of the area. If failure strengths 
are used, then the aforementioned parameters become 
variables because of variations in weight, and the more 
simple problems such as, for instance, pure tension or 
compression under unit probability of loading become 
considerably more complicated. However, it will be 
seen that in the general case, which is treated in this 
paper, where the probability of loading is taken into 
account and where the failure stresses are dependent 
on the area, the standard deviations and mean stresses 
turn out to be unknown variables. 

For the sake of simplicity let the applied and failure 
stress distributions be normal, such that 


(fm’ —fm)?/2afm? (1) 


-— r .— 
mT m ) _ (1 270 7, )é 


ai / ‘ ,— (Fm! — Fm)?/20Fm? - 
el Tas ) = (1 V/2ror,)é ” (2) 


The assumption of normal distributions is accompanied 
by certain limitations. For large values of the standard 
deviation (15 to 20 per cent of the mean and higher), 
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MINIMUM 


the normal distribution curves are such that they have 
relatively large values in the negative region of the 
abscissa. This would imply that for low probabilities 
of failure it would be possible to obtain, say, tension 
failures due to compressive loads. Such a phenomenon 
has not been observed, indicating that the present 
analysis must be limited to standard deviations of 
allowable stresses of less than about 15 per cent of the 
mean. This is not an unreasonable restriction since it 
is highly unlikely that a material would be selected 
with the wide scatter of allowable stresses correspond- 
ing to these high values of standard deviations. For 
higher values of standard deviations, distribution 
curves which are zero for negative values of the variable 
must be chosen. 

The use of Gaussian distributions, however, does not 
detract from what is to be derived, i.e., the laws of dis- 
tribution of failure probabilities of the various struc- 
tural elements such that the total structural weight is a 
minimum and the probability of failure of the entire 
structure is equal to a preassigned value. Since this 
is the first time such an analysis has been carried out, 
Gaussian distributions have been chosen to simplify 
the computations. The such distributions 
affects only the numerical results and does not influence 
the rational design procedures to be formulated. 


Let the difference between the failure and applied 


use of 


stresses be g, 1.€., 
San = Fa’ = ee (3) 
By Egs. (1) and (2) it can be shown that g,,,’ will also 


be normally distributed with a mean value g,, and 
standard deviation ¢,, such that 


£m = ry boise? ie (4) 
Om” = OFm” + Fy" (5) 


Failure will occur whenever f/f,’ > /,,’, so that the 
probability of failure of any element is found by inte- 
grating over that region of g where it is negative. 


Therefore 
0 1 g—2 2 
P,, = (1/V 2a on) f e ? ( or ) dg (6) 
Let x = (g — BZm)/On (7) 
Gn - £m Tm (S) 


Then Eq. (6) may be written in the form 


Ge 
P, = 1/2 =— (1 V2) { e~*'* dx (9) 
0 


Thus it can be seen that the probability of failure of the 
element is a function of one variable only, G,,, which 
now must be related to the weight and probability of 
failure of the entire structure. The latter is given by 


the expression 


P,=1-—(1—P)(1 — Pe)---(1 — P,) (10) 


The probability of failure of the entire structure, P,, 
is prescribed by specification or other choice, and it is 
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generally small, say, of the order of 1 in 1,000. The 
probabilities of failure of the individual elements mak- 
ing up the entire structure must obviously be less than 
P,. Consequently, one can neglect all nonlinear terms 
in Eq. (10) and obtain a good approximation from 


P,= 2 Pr 
m=1 


(11) 


Substitution of Eq. (9) into (11) leads to 


Gm 
> e~*"? dx = +/2 (° — P.) (12) 
m=1J0 2 
It will be assumed that the length of each element is 
fixed in the design of the structure and that only the 
individual areas can be varied to achieve minimum 
weight. It will be shown subsequently that both the 
mean failure and mean applied stresses and standard 
deviations of each element can be expressed as func- 
tions of the corresponding areas A,,. Since the prob- 
abilities of failure are dependent on these quantities 
through Eqs. (4), (8), and (9), one. may express all 
variables in terms of the respective areas. It then 
follows that Eq. (11) may be written in the form 
P, = > P»(Am) (13) 
m=1 

The weight of the entire structure, IV’, is the sum of 

the weights of the individual elements: 


W = D Wa = & nlmAm 
1 


m=1 m 


(14) 


The purpose of this analysis is to minimize the total 
weight IV. subject to the restraint of Eq. (13). This is 
best accomplished by the use of Lagrange’s method of 
undetermined multipliers,* which results in a set of 7 
equations of the type 


OW/0A m + w(OP;/OA») = 0 (15) 


where yu is the Lagrangean multiplier. 

The solution of the system of Eqs. (13) and (15) is 
given in detail in Appendix B. Basically it consists 
of the elimination of the Lagrangean multiplier yu 
between pairs of Eqs. (15). However, before the solu- 
tion is presented certain facts pertinent to the deriva- 
tion must be discussed. 

In Appendix A it is shown that for a large class of 
loadings the mean stresses can be represented by the 
relations 


i = nas Am A m (16) 


Fm = KmAn*” (17) 


where X,, is the mean applied load of the mth structural 
element. The constants K,, and N,, are tabulated in 
Appendix A and are functions of material properties 
and the type of loading. The results to follow are 
derived in terms of these constants. 

For the purposes of this analysis it is convenient to 
introduce the parameters 
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Crm TF m Fe (18) 
Crm = CO fm) j= ( 19) 


It will be assumed that for any given region of F,, and 
fm, there correspond constant values of Cy, and Cy, 
which are known from prior tests determining the 
statistical failure stress distributions, and from tra- 
jectory analyses giving load distributions. 

This approximation is equivalent to replacing the 
presumably smooth functions C,,,(/,,) and Cy,(fm) by 
an arbitrary number of step functions. For N,, = 0, 
(tension, compression, etc.) Cy,,(/,,) is indeed constant 
for a given material. The test data of Harris et al.’ 
indicates that C,r,, is nearly constant over a wide range 
of F,, values for the buckling of thin-walled cylinders 
(N,, = 1). 

For applied stresses of the type defined by Eq. (16), 
one can show’ that C,,, is related to the standard devi- 
ation of the applied load (o,,,) and of the area (¢,4,,) 
by the approximate relation 


a ( OA m A 2 ( 20) 


‘ ° ” 
Cm" = (Orn, Xm)” 


The ratio o4,,/A» is a function of the manufacturing 
processes and has been found to be a constant. No 
data seem to be available for the statistical distribution 
of loads in a given flight; however, it will be assumed 
that either o,,,/A» is a constant or that it can be ap- 
proximated by step functions over various ranges of 
Am. Hence, the assumed behavior of C,,, and Cr,, are 
similar in character, and they are treated as parameters 
in the solution to follow. 

The details of the solution of Eqs. (15) are presented 
in Appendix B and the results are symbolized by Eqs. 
(21), (22), and (23): 


Zz Lip N; + 1) /r, 1/(Ni+1) K\! (Nj+1) 
oe SGP V i j = 
Z; : Lip(Ni + 1) (i) (¢ ) 


(t Cr2G?)'*”) (Ni+1) (1 + GB) (Vi+1) 


FE Vit, : ——— x 
i-(A7 "CC hla” 
(2: + cue aia 21) 
3 —=f~s0/@ (Z 
B, + Cr?7G,/ Bye~°"” 


1/(Nm+1) Y 1/(Nm+1) 
W I ¢ ) ( l + Cuban ) vs 
m = Pm sm - Pe ’ 97. 9 
K m 1-—C Fm °G m 


(22 


WN, + 1) 
W(N, + 1) 


( = ret L+ Bs) = Be" 
_ Cr/*G;" 1 + G,B,/\B, + Cr’G; Bye o* 


(23) 


where Ba? = Cm? + Crm? — (Chm CrnGm)* (24) 


It should be noted that these results appear as functions 
of properties of pairs of structural elements. 
tion (21) relates the difference between the mean failure 
and applied stresses, G,,, of any two elements, to 
known parameters of length, density, and type of load- 


Equa- 


1969 
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ing. The quantities G,, are, of course, related to the 
probability of failure P,, of each element through Eq. 
(9). One may therefore think of Eq. (21) as a relation- 


ship between probabilities of failure of pairs of elements. 
' 


Equation (23), on the other hand, relates the G,, and, 
consequently, the probabilities of failure of pairs of 
structural elements, to their respective structural 
weights. Equation (22) gives the relationship between 
structural weight and probability of failure of each 
structural element. A more detailed discussion of these 
results will be given in a subsequent secticn. 


The Margin of Safety 


It is customary in present-day analyses, where proba- 
bility concepts are not used, to compare design stresses 
with allowable stresses in terms of the margin of safety. 
Such a procedure seems to inspire some degree of con- 
fidence. However, the mere fact that one may claim 
that a certain structural element has a given failure 
stress and a preassigned design stress does not make it 
so in reality. For instance, in one thousand aircraft a 
structural element will have a range of statistically 
distributed failure stresses, and each aircraft will be 
exposed to loading conditions in which the occurrence 
of this one particular design stress is only a statistical 
probability. 

Obviously, if the probability of failure of this struc- 
tural element is, say, one in a thousand, then failure 
will occur with this probability regardless of whether 
a positive margin of safety is achieved in the analysis. 
The margin of safety is based on a single combination 
of design and allowable stresses. On the other hand, 
the thousand structural elements under discussion here 
may have a thousand different such combinations. 
On the basis of the one particular combination used in 
the margin of safety equation, the computation may 
show a positive or negative answer, but the numerical 
value of the margin of safety, in this instance, bears 
little relation to the failure of each of the thousand 
elements. 

However, in order to tie together the present-day 
method of analysis and the concepts of structural 
reliability used in this paper, relations are derived here 
between probability of failure and the margin of safety. 
At the outset it must be pointed out, in all fairness, that 
the authors have used a certain prescription for select- 
ing a pair of design and allowable stresses to compute 
numerical values of the margin of safety. While it ts 
true that the pair of stresses used here in the computa- 
tions is the currently accepted combination, others 
may wish to lay claim to different pairs of such stresses 
and arrive at different numerical values of the margin 
of safety for the same probability of failure. The ex- 
pressions derived herein are general enough to permit 
such selections, but one must keep in mind that a 
particular numerical value of the margin of safety is 
inherent in the selection of a particular pair of stresses. 

The margin of safety of the mth element is defined 
in the usual manner: 
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(M.S.)m = = — (25) 
where the design stress can be related to the mean 
applied stress by 


Tar - lm 2 Thm. (26) 


and the ‘‘allowable stress’ F,,” can be related to the 


actual mean failure stress by 


; ' - 
F.” = F, — ore” (27) 


" 


It is customary to take F,,” as the minimum guar- 
anteed value of the failure stress; however, for the 
purposes of this analysis any preselected value will do 
The design stress is usually selected from 


equally well. 
the worst loading condition, without regard to the 
mean applied stress. A factor of safety may be intro- 
duced into Eq. (25) at will by the proper adjustment 
of the value of a,,,”. 

Substitution of Eqs. (18), (19), (26) and (27) into 
Eq. (25) and division by ¢,,, yields 


° ] s o fm l es 
(M.S.), = : — X Pus : ; + X Fs — | 
( Fm TFm ie Im 


(28) 

where Xp, Or,” / Orn (29) 
X fn = Fjm” | F tm (30) 

1 For a given probability of failure Cy,, and Cy, are 


known, and o,,,” and a,,,” may be chosen by the de- 


signer. Therefore, only the ratio o,,/or,, must be 
determined in terms of known values. 
By use of Eq. (19), Eq. (5) can be expressed in the 


form 





(OF m/F jm)? = (Om/Cimfm)? — 1 (31) 


and substitution of Eq. (89) yields 
O t,/ O Fm Cal l1-— Cetin) Cr,,| | + ByGm) (32) 


From Eq. (28) and (32) one obtains 
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Fic. 3. Normal distribution vs. probability of failure 


. i = oe Oe l =< 'F Nn X Fp 
(M.S.)m = nn —— —1 (33) 
1 — Cr,%Gn'/\1 + Cy Xy, 


Consider next the case where F,,” is taken as the 
The latter 
usually expressed as the lowest stress expected in J, 


minimum guaranteed value. value is 


percent of the material. As in the previous section, 
a normal distribution is assumed, and thus 


Pym 


Let y= (F- (35) 
then with the help of Eqs. (27) and (29) one obtains 
from (34) 


XF 
Vn = (1/2) + (1/V 2r) [ e** dy (36) 
7 


which gives the necessary relation between the known 
Y,, and the corresponding X ,,,. 

The relationship between the margin of safety and the 
probability of failure is shown in Fig. 11 and will be 
discussed subsequently. 


Weight Saving 


In order to evaluate the benefits of the analysis con- 
tained in this paper over conventional analyses based 
on margin of safety considerations, a comparison is 
made between the weights of structures computed 
For the sake of simplicity 
we will consider a two-element structure. Let I") and 
I. be the weights of the respective elements as com- 


from each type of analysis. 


puted on the basis of margin of safety, and let P; and 
The 


comparison will be made under the following stipulated 


P, be the corresponding probabilities of failure. 


conditions: 
(1) In the conventional structure, the margins of 
safety of all elements are equal. 
2) The probability of failure of the conventional 
structure is equal to that of the optimized structure. 
For the two-element structure the weight saved, 
AW, is given by 
AW/W = 


(, + WT. — (Wi, + W2))/(Wi + We) (37) 
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Rearranging Eq. (37) in order to substitute in it terms 
of previously determined quantities, we obtain 


AW/W = 
W(/W 2) + 1 /W2[(W4/W2) + 1] — 1 (38) 


Use of Eq. (22) yields 


where 


— | ( 1 + GB, N° iia ca) (Ni +1) mm 
os 1— Cr2G? 4 Ge. ( 

C ( 1 + GE. \(? = CnlG\T (N+) - 
= @ ( 
2 l1— Cr,°G? 14 oh 


Substitution of Eqs. (22) and (89) into (38) gives 


AW/W = [C\(W/W2) + GC) (Wi/W2 + 1] -— 1 
(42) 


The ratio V/V. can be determined from Eq. (23) for 
a given set of values of P,, Cr,,, C,,, and NV. 

The barred quantities in C,; and C, are determined 
from the two conditions prescribed at the beginning 
of this section. Use of Eq. (33) gives the following 


relation for condition (1) 
( 1+ BG, (; — cnet) 7 
l _— Cr2*Gi" ] + Cy, X y, 


( 1+ BG, )( — cele) - 
= (+9 
1 — Cp,?Go?/\1 + CyX 4 


and Eq. (12) determines condition (2) to be 


6, G / 
J e7*/? dx + f ee"? dx = V2r(1—P,) (44) 
‘ 0 


Eqs. (43) and (44) are sufficient to determine G; and 
G, in terms of other known quantities. 

The results of the numerical computations pertaining 
to weight savings are shown in Figs. 12, 13, and 14 
and will be discussed in the next section. 
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Discussion 


Numerical computations for a relatively large num- 
ber of values of the several variables were perforined on 
an IBM 704 digital computer at Hughes Aircraft 
Company. However, because of space limitations 
only the most significant results will be discussed here. 

As was shown previously in the Analysis, the solution 
of the minimum weight problem consists of the simul- 
taneous solution of Eqs. (12) and (21) for the un- 
knowns G,,. To facilitate this procedure, Figs. 2 and 
4+ have been prepared to be used directly in conjunction 
with Eqs. (21) and also with Eqs. (22) and (23) for the 
calculation of the structural weights of the various 
elements. 

From the calculations that were made. it was ob- 
served that Eqs. (21) and (23) could empirically be 
reduced to simpler forms which, as will be seen, will 
greatly reduce the complexity of the many-element 
problem. The numerical results, such as those shown 
in Figs. 4 and 5, indicate that Eqs. (21) can be written 


as 
Pi/ Pj = a4jQij"" 45) 
where 


Q a (Z; Z (Coe? + Cora 7 x 


(Cy? + Crj7]-"? = (2,/2;)Cj* (46) 


and, similarly, that Eq. (23) can be reduced to 
P,/P, = 6.2." (47) 


with 0, = Ciy* WAN; + 1)/Wi(N; + 1) (48) 


Sij 

Eqs. (45) and (47) are basically simpler equivalent 
approximate expressions of Eqs. (21) and (23), and as 
such may be preferred by the designer. 

Figures 4 and 5 were obtained from Eqs. (21) and 
(23) for a representative set of values of P; + P,, 
Crim» Crm With N; = N; = 0. The constants a;;, 5 
a;;, and £;;, also functions of the above parameters, 
were empirically determined by the method of least 
squares from the computations shown in Figs. 4 and 5 
and are presented in Figs. 7-10. The radical C,,* in 
Eqs. (46) and (48) was introduced as a scaling factor, 
because it was noted that its use brought the curves 
in Figs. 4-10 closer together. 

Figures 4-10 represent cases where the largest 
spread of results was encountered. The largest differ- 
ences between the approximate results of Eq. (45) and 
the exact results obtained from Eq. (21) are 3.39 per 
cent for the Q curves and 2.76 per cent for the 2 curves. 

It can be seen from Eqs. (21) and (46) that for any 
structure Q;,; is always known since it depends solely 
on the length, density and type of loading of pairs of 
structural elements. In some cases it may be desirable 
to know Q,,; directly without computing P;/P, before- 
hand. This ratio can be eliminated between Eqs. 
(45) and (47) which results in 
Qey = (Vy ae/b,) 1% (49) 


**ij 


A plot of Eq. (49) is presented in Fig. 6. 
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Certain important conclusions can be drawn from 
Eqs. (47) and (48). It can be seen that whenever 


bij [(Cij*(N, + 1)/(Ni + I)]8" < 1 (50) 


then W,/W; > 1 for P;/P;> 1, or the heavier element 
has the higher probability of failure for minimum over- 
all structural weight. The computations show that 
even when the inequality (50) is reversed, that as the 
component becomes heavier it should also be given a 
higher probability of failure for minimum overall struc- 
tural weight (see also Appendix C). 

The aforementioned empirical results may be used 
to good advantage in the solution of the many-element 
structure. The following procedure may be employed. 
Equation (11) is rewritten as 


P, = P(i + > Pn p,) (51) 


m=2 


Then, solving for P; and using Eq. (45), 


; 1 
P, = pli + am Oms"™ | (52) 
m=2 

Unfortunately, since each a,,; and a,,,; depend on the 
sum P,,,, = Pm + Pi, an iterative process must be used 
to estimate pairs of @,; and a, values. For the first 
cycle one assumes values of /,,,, and obtains the a,,, and 
a, from Figs. 7 and 8. The Q,,1 are known, of course, 
since Eqs. (21) and (46) define them exactly. Then P; 
can be found from Eq. (52) and the values of P,,, can be 
found from Eq. (45). If the P,,,, values now obtained 
differ materially from those originally assumed, then 
one may repeat the computation of P; by using new 
Qmni and @m, values corresponding to the P,,,, just 
found. This process may be repeated until satisfactory 
convergence is reached. 

From the last value of P; one may obtain G; from 
Eq. (9) or Fig. 2. Then one uses Eq. (22) to compute 
W,. The weights of the other structural elements are 
found from Eqs. (47) and (48) and are given by 


Wm = WiCm*((Nim = 1)/(My = 1)] t= bm P| ' bi 
(53) 


SCIENCES SEPTEMBER, 1960 
The 4,,, and 8,,; values must, of course, correspond to 
the previously determined P,,; values. 

Figure 11 shows the variation of the margin of safety 
with probability of failure for representative values 
of X,, Y, Cr and C; as computed from Eg. (33). It 
has been pointed out already that the numerical value 
of the margin of safety is misleading since even a posi- 
tive margin of safety does not give assurance of a com- 
pletely safe structure. The trend of margin of safety 
values is as expected, 1.e., the lower the probability of 
failure, the higher margin of safety. The numerical 
values of the margin of safety in this paper are based on 
mean applied loads and the design safety factor of 1.5 
has been omitted. 

It should be noted from Eq. (33) and Fig. 11 that 
the mere selection of design values of X »,, and X s,, does 
not ensure a unique value of the margin of safety. The 
margin of safety depends on the probability of failure 
through the latter’s relationship to G,, [Eq. (9)|. 
However, the probability of failure is independent of 
the relative dispersions Cr, and Cy. Consequently, 
for any Xp,,, X y, and P,,, various values of the M.S. can 
be obtained depending on what statistical material 
properties are involved (Cr,,) and on what the flight 
conditions are (Cy,,). 

It is of interest to note that the proportioning of 
probabilities of failure to the various structural ele- 
ments as shown above is equivalent to the distribution 
of margins of safety to these same elements. Use of 
Egs. (21), (22), (33), (35), and (48) yields 


0, - (: _ Cook | ] Vit] , 
S ti 4+ ¢.%27 (MS +1 
| Gs) M.S.) +n] ah 54) 
( VE wD.) 7 () 
1 — CrjXr, 


Then from Eqs. (45), (47) and (54) one obtains 


<: lias ter f 1 \™ 
(By aves (1) x 
l — Cr:X ri | 31/(Ni+1 
( + CX pi (M.S.); + | x 
1+ C,Xy, W/(Ni+D 
(| —c x) }(M.S.); + u | (53) 
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which gives the equivalence between the margin of 
safety and probability of failure ratios. 

Finally, substitution of Eq. (55) into Eqs. (45) and 
(47) results in 


[M.S.),; + 1)/%* [(MLS.), + IT VOtY = 
( + CX gi ) “( _ serie (Ni+1) . 
lL — CriXr; L + CyX 4; 


b,, 1/Bij 
( 0," aij — 1/Bij (56) 
Qi; 
and 


[((M.S.), oo 1]! Vit) 1(M.S.); + we (Ket) 
( + CyiX 7 ) ~G 8 cen (Ni+)) " 
aa CriX Fr: l + CyjX yj 


l/a 
b; sg 
j 1/aij — 1/Bi == 
( Qi, » kee) 
Qi; 


A comparison of Eqs. (56) and (57) with Eqs. (45) and 
(47) shows that the latter (involving the probabilities of 
failure) are considerably simpler and more direct, 
much easier to compute, and easier to visualize. In 
addition, the numerical values of the margin of safety 
ratios depend upon the arbitrary selection of Xr and 
X, values, while the probability of failure ratios have 
unique values for any Qj). 

Figures 12 and 13 show the weight saving, AW/W, as 
a function of the weight ratio W,/W,;. These curves 
were computed by using Eq. (42) for representative 
values of Cr», and Cym, for one set of X,and Y values, 
for N, = N,; = 0 and for P,,;, = P; + P; equal to 
10-* and 10-8. As indicated in the figures, there is 
for any combination of Cr», and Cym a zero point, i.e., 
a value of W,/W, which is equal to W,/W,, where the 
weight saving is zero. The weight saving increases 
to asymptotic values, which are independent of W;/W,, 
as W,/W, increases or decreases from the zero point. 

As seen in Fig. 12, the range of weight savings is 0 to 
19 per cent for P,,;, = 10~*. For the same values of 
X,and Y, the range of AW/W increases with decreasing 
P,,,, reaching a value of 74 per cent for P,,;; = 107%, 
as illustrated in Fig. 14. 

The average weight saving in a multiple-element 
structure will generally be less than the maximum. 
For the simple model, with a Gaussian distribution and 
a two-element structure, the weight saving may vary 


from 0 to 74 per cent, depending upon the values of the 
various parameters NV», X jm, Ynys Cm, Crm and Pj; 
The weight saving for a multiple-element structure 
will be between the minimum and maximum. In order 
to estimate what the average value might be, the fol- 
lowing technique was used. The average value of 
AW/W for 10,000 two-element comparisons for various 
values of the above parameters ranged from 2 per cent 
to 8 per cent for P,,, equal to 10~* and 10~* respectively. 
This averaging process, shown in Fig. 14, may give 
an indication of the average weight saving of the entire 
multiple-element structure. 

Figure 15 is a qualitative illustration of the weight- 
saving mechanism; case E of Fig. 12 is used as an ex- 
ample. The solid curves represent a state based on 
conventional analysis where the total probability of 
failure is 10~* and the margins of safety of the two 
elements are equal. However, P,/P; is very large, i.e., 
maximum weight savings occur whenever either P; or 
P, is nearly equal to the sum P,; + P, (in this case 
P,= P..). The dashed curves represent the minimum 
weight for P; and P; determined from Eqs. (43) and (44) 
and (9), their sum being again equal to 10 Large 
weight savings can be achieved by a small decrease in the 
applied stress in the element 7, which has a probability 
of failure nearly equal to P,,,, and by a large increase in 
the stress of element 7. This results in a small decrease 
in P; and an increase in P;. Thus, even if W; = W,, 
the net weight saving can be roughly approximated by 
an increase in applied stress or a decrease in area of the 
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element having the smaller probability of failure, i.e., 
element 7 in this case. 

All the results presented in this paper are based on 
Gaussian distributions, though admittedly, actual dis- 
tributions are nonGaussian. However, the choice of 
distributions affects only the numerical results and not 
the general pattern of the problem. Research is cur- 
rently in progress to extend the analyses and design 
procedures to nonGaussian distributions. 


Appendix A 


Derivation of General Relations Between 
Failure and Applied Stresses and Weights 


In this section Eqs. (16) and (17) will be derived for 
various types of loadings. The loadings treated here 
are, of course, not all-inclusive since Eqs. (16) and 
(17) do not apply for inelastic structures and for com- 
bined loadings. In particular, the statistical distribu- 
tions of allowable stresses in combined loading, such as 
peel stresses in bonded joints,’ do not have mean failure 
stresses representable by Eq. (17). Because of the 
inherent nonlinearities associated with the combined 
load problem, a somewhat different approach must be 
used which may be found in reference 10. 


Simple Tension or Compression and Tension 
or Compression Creep 

For this type of loading X,, represents the load in the 
element. The mean failure stress is independent of 
the area and is equal to the mean ultimate stress. 
These results can also be applied to creep failures: for 
preassigned lifetime and creep strain, there corresponds 
one mean stress which is independent of area. 


Elastic Buckling of Circular Cylinders 

Here X,, is defined as above; however, the mean fail- 
ure stress is equal to the mean critical stress which is a 
function of the area, since 


Fy = AnEgte/Ru (58) 
For a circular cylinder the cross-sectional area is 
An = 24Rntm (59) 
and substitution of Eq. (59) into (58) yields 
F,, = (@,,/27R,*)A» (60) 


Generally, the radius R,, is fixed because of other de- 
sign considerations, for instance, the size of the fuselage, 
and only the thickness ¢,, may be optimized. Then 
all the terms except A,, on the right side of Eq. (60) 
are constant and equal to K,,. 


Creep Buckling of Circular Cylinders 


Bozajian and Tsao! have shown that the critical 
or failure stress due to creep buckling of circular 
cylinders can be predicted from 


(61) 


»' . / 
Pins = BB albert ai Rn 


AEROSPACE 


1960 


SCIENCES—SEPTEMBER, 
where the failure probability factor 8,, is a function of 
the probability of failure P,,,, or 


Bm sa BoP.) = r* Pas (62) 


The secant modulus, £,,,, is a function of the maximum 
time, 7 mar, at buckling and the applied load X,,’ 


Bo = Bos + tea 4 (63) 


m 


where 


Teas k 
Cm aaa 1 Emm [\m’(Tmaz)] | | f Dive (r) |" . ar (64) 
0 


The values of the constants b,,, , and k are given in 
reference 11. From Eqs. (59), (61), (62) and (63), one 


obtains 
F,, = 4,8 Aw) (2eh,*(1 + teAg! yi (65) 


An examination of Table 1 of reference 11 shows that 
for elevated temperatures the inelastic strains are 
considerably larger than the elastic ones. In terms 
of the nomenclature of this paper, this means that the 
term ¢,,A,,'~—" is much larger than unity and Eq. (65) 
can be approximated by 


Fry = dmEnA m"/(247Rm¥m) (66) 
If one treats ¢,, as a parameter that depends, as was 


seen previously, upon the loading history and failure 
time, then Eq. (66) is also of the form of Eq. (17). 


Elastic Buckling of Plates 


In this instance X,, is defined as previously and the 
mean failure stress is again equal to the mean critical 
stress, which is now 


P= oa /13l = v7) / On)" (67) 
The loaded area of the plate is 
A. = 6S. (68) 
which reduces Eq. (67) to 
Fy = knEint*A»?/(126,4(1 — v»,")] (69) 
If one assumes that the size of the plate is fixed by other 


design considerations such as stringer spacing, then 
both &,, and b,, are constants. 


Bending or Torsion 

In the case of bending or torsion of circular tubes, 
the mean allowable stress is the mean modulus of rup- 
ture in bending or torsion which can be approximated 
by step functions and considered as a set of different 
constants for small variations in area. 

The mean applied stress in bending is 


Pag = DM glad Ti (70) 
and in torsion 


Ne = 7 Pr m J, m ( 7 ] ) 


where /,, and 7,, are the mean bending moment and 
torque respectively. The moments of inertia are 
related to the areas by the radii of gyration, such that 
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In =A (73) 
which modifies Eqs. (70) and (71) to read 


Sm = MaCm/(Qm7Am) (74) 


and 


lc = Tula Cte ad (75) 


For a fixed outer radius c,,, the radius of gyration may 
be assumed to be constant and, consequently, one opti- 
mizes only the thickness of the tube. Then the first 
group of terms on the right-hand side of Eqs. (74) and 
79) may be designated as \,,, since they vary with 
applied moments and torques only. 

Table 1 summarizes the results obtained in this sec- 


tion. 


Appendix B 


Derivation of Equations (21), (22) and (23) 


In the Analysis it was stated without proof that the 
solution to the problem was given by Egs. (21), (22), 
and (25). These relations can be derived in the fol- 
lowing manner. 

First one substitutes Eqs. (13) and (14) into Eq. (15), 
which results in 


or. . a 
= 0 (eka 22S. ws ct) tre 


ie ~ 
, ” ad. 


Equations (13) and (76) constitute a system of r + | 

equations in unknowns (u and r number of A,,). The 

multiplier « can be eliminated between pairs of equa- 

tions of the system (76), which are characterized by 
m = 1andm = j, such that 

pil.i/ pj; = (OP;/0A;)/(OP;/0A;) (77) 

Finally, use of Eq. (9) in (77) yields 

pil; (e~ °*) OG,/0A; 

pyL, (e~ %"*) OG,/0A, 


(78) 


There remains now only the determination of the 
derivatives O0G,,/0A,. From Eqs. (4) and (8) one 
finds that 


G, 1 ES Wn _ on | on 
uM” oie. & ° "1 


From Eqs. (16) and (17) one obtains by differentiation 


and proper substitution 


Of m OA eae —fim? Mini (SO) 


dF ,,/dA m = fmE mN m/ Ym (S1) 
Using Eqs. (18) and (19), we can now write Eq. (5) 
in the form 

Om? = Crm? Fm? + Cim?fm? (S82) 
Differentiation of Eq. (S2) and use of Eqs. (SO), (81), and 


82) yields 
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TABLE 1 
Summary of Relations for the Equations 
fm = dAm/Aw and F, = Kn Am** 
Loading Am K, N 
Tension or Mean applied Mean ultimate stress 0 
compression load 
Tension or Mean applied Mean stress correspond 0 


compression load ing to preassigned 
creep time and creep 
deformation 

Elastic buckling Meanapplied a,,E»/27rR,? l 
of circular load 
cylinders 

Creep buckling Meanapplied  d»,Egm/27rRmt, n 
of circular load 
cylinders 

Elastic buckling Meanapplied k,E»,7?/12(1 — yp?)b»!4 2 
of plates load 

Bending of M m€m/2m? Mean bending modu 0 


round tubes lus of rupture 
Mean torsional modu- 0 


lus of rupture 


Torsion of round = 7n€m/2nm? 
tubes 


dom dA m ( Cala Am) x 


[Nm — (Nm + 1)(O4m/Om)*] (83) 


Substitution of Eqs. (80), (81), and (83) into (79) now 


gives 


OGm/0A m = (Nm + 1)(fm2/Omdm) X 


(1 + GinCgm?fm/Om) (84) 


One proceeds next to the evaluation of f,, and o,, in 
terms of known quantities. Substituting Eqs. (16) 


and (17) into (82) and (&) results in 


— : no ee fre Sets — ‘ 
(Om/fm)? = [Cra?K m7\m?/f mn |] + Cy,” (85) 


Vmn+1 » No 


(Om/fm) = (KmAmi" — fm \/fm"** Gm) (86) 


Squaring Eq. (86) and eliminating (¢,,/f,)* between 
this new relation and Eq. (85) gives 
2 > 


» = finde oo SNe oe 
| = PR (1 — G,,*¢ Fn) = 


fa Set (20K dm") +1 — Cy2?Gn? = 0 (87) 


This is a quadratic equation in /,, and has the 


usual solution 


fim NX" = KadmX™ (1 — Crm?Gm?)/(1 + GmBm) (88) 


where B,, is given by Eq. (24). 
From Eqs. (86) and (88) one obtains 


(Om/fm) (Bm + Crm?Gm)/(1 — Gm?Crn”) (89) 


The positive sign in the denominator of (88) and in 
front of B,, in (89) can be justified in the following 
manner. Consider the case where C,,, Ce... Ac 
cording to Descartes’ rule of signs,'* Eq. (87) will have 
positive real roots only if 


GCs.* = 6'C,* < 1 (90) 
Then one obtains from Eqs. (24) and (89) 
B,* = Cr,,°(2 — Cr,,"Gu*) > Cra’ (91) 
and from (90) 


CoG < Co (92) 
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Therefore, from Eqs. (91) and (92) one concludes 
Bm Cain (93) 


and the positive sign must be chosen if Eq. (89) is to 
remain positive. A similar proof can be offered for 
the cases where Cr,, # Cy,,. 

From Eqs. (88) and (89) one obtains 


4 Xr 1/(Nm+1) 
( f,.2 a) = = m x 
; (, + s_ | 


(Nm+2)/(Nm+1) 


(1 — Cpr,,°G,,”) 


= (94) 
Bes + ( a, oe 


and substitution of Eqs. (89) and (94) yields 


0G,,/0A m = Bu(Nm + 1)(Km/Am) Ot” & 


’ val (Nm+2)/(Nm+1) 
(1 — ¢ Fin ig) 


Ee (14+G,B,)°" "rt" (95) 
(5. 4 Cetra) 


which completes the evaluation of all derivatives in 
Eq. (78). 

Using the results of Eq. (95), one may now express 
Eq. (78) in the form of Eq. (21). Admittedly, Eq. (21) 
looks complex but it must be remembered that it is 
only an algebraic equation in the unknowns G; and G;. 
The solution of the problem for the structure consisting 
of r elements is the simultaneous solution of r — | 
equations of the type (21) together with Eq. (12) for 
the unknown G,,. It should be noted that the left- 
hand side of Eq. (21) contains only known quantities 
of structural element length, density, and type of 
loading. The quantities Cy and Cy, are parameters 
which are known from material and loading statistical 
distributions. This system of equations can be 
solved numerically as was discussed previously. 

Eq. (22) is obtained from Eqs. (14), (16), and (88), 
and substitution of Eq. (22) into (21) yields Eq. (23). 


Appendix C 


The Relationship Between Weight and 
Probability of Failure of Structural Elements 


It is of interest to note that it is possible to set down 
certain rules for the distribution of probabilities of 
failure among structural elements such that the overall 
weight of the structure is a minimum. 

For the sake of simplicity, consider the case of unit 
load probability, C;, = 0, for which Eqs. (21) and 
(23) are considerably simpler: 


Z. f #NEED/ND IL pp #1/(NIFD IL 9 ik Fit)? 
i Ji J ; 
Z, a | Seevieeeeren E wma | aa | 
JJ 1 


(96) 
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f;* — F{*)2 


WAN, +1) _ fit e7 


W(N, +1) f,% eo 9*-F (97) 

with 
fm* = fn/(W2 oFm) (oa) 
Fn* = Fu/(/2 ov.) = 1/(/2 Cr.) i 


Consider Eq. (97). If WAN; + 1) > W,(N; + 1), 


then it follows that 


f2_  ¢g GP-FIP gt”) 

ci” ae Se (100) 
dj e @;(fi*) 

Now, if ¢)/¢; > 1, then from Eq. (97), f;*/f;* > 1. 
However, because of the relationship between the ¢’s 
and f*’s [Eq. (100)], these two inequalities are incon- 
sistent. On the other hand, if ¢;/¢; < 1, then /;*/f;* > 
1 and the inequalities are compatible. 

Since the coordinates of the normal distribution 
curves are related to the probability of failure of the 
element, one may conclude from the above that for 
N,; 2 WN, the heavier element will have the higher 
probability of failure. 
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Terminal Guidance System for 
Satellite Rendezvous 


W. H. CLOHESSY* anp R. S. WILTSHIRE** 
Lhe Martin Company 


Summary 


This paper assumes a requirement for an unmanned multiunit 
satellite to be assembled in orbit. The requirement to be met is 
to bring the satellites together so that they do not collide but 
actually rendezvous. The equations of motion of the rendezvous 
satellite in a relative coordinate system are derived and used to 
compute a final injection velocity which would effect collision 
after a time r. The velocity is corrected periodically by a com- 
mand guidance system and just before impact retrothrust is ap- 
A terminal infrared homing system is required to actually 
The 
first satellite placed in orbit is the ‘‘control satellite’? and controls 
all the satellites to be assembled and contains the cc mputer, 


plied. 
accomplish physical contact and joining of the satellites 


command guidance equipment, precision orientation equipment, 
and other features necessary to effect rendezvous. The succeed. 
ing satellites contain a propulsion system, a rough attitude con- 
trol system, and a command receiver plus whatever scientific 
equipment they carry to perform their basic mission 

This paper presents the following: 


(1) A description of how the “control center’s’’ orbit is deter- 


mined and made available to the ‘‘control center’? computer. 


2) The attitude sensing and control system of the 


(2 


center” and of the ‘“‘rendezvous”’ satellite. 


“control 
(3) The measurement system by which the ‘control center’ 
determines relative position and velocity of the ‘‘rendezvous”’ 
satellite. 

(4) The equations which must be solved to generate the 
rendezvous commands and the computer used to mechanize these 
equations 

(5) The command link from 
satellite and the receiver and internal control system on the 


“control center’’ to ‘‘rendezvous”’ 


“rendezvous” satellite. 

(6) The homing guidance system to effect closure of the last 
few hundred feet. 

(7) A description of the dynamics of the rendezvous itself. 


Introduction 


i ies ‘“‘RENDEZVOUS”’ PHASE of the process of building 
a space station or of assembling, in orbit, a big 
satellite out of little ones is the subject of this paper. 
More particularly the dynamics of ‘‘rendezvous’” and 
the instrumentation necessary to bring about ‘‘rendez- 
vous” are discussed. The orbit altitude is not severely 
restricted but has been chosen such that an integral 
number of revolutions occurs in 24 hours, almost. Ac- 
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tually because of the precession of the nodes of the orbit 
the altitude is shifted so that the orbit plane will contain 
the launch sight, exactly, every 24 hours. The succeed- 
ing satellites or rendezvous satellites could be fired at 
one day intervals, burned out at some positive angle to 
the local horizontal transfer orbit, and allowed to coast 
up to apogee at the altitude of the control satellite. 
This manuever is timed so that apogee is reached at a 
position and time near the arrival of the control satel- 
lite and perhaps slightly ahead of it. How actual ren- 
dezvous is effected will be shown later. 


Tracking and Stabilization 


The initial satellite (the ‘control center’’ or ‘“‘control 
satellite’’) is assumed to be in orbit, in a circular orbit. 
It is tracked from the ground by means of a radar. 
For this purpose we have chosen an FPS-16 (RCA) 
radar whose characteristics are described in Fig. 1. 

The required signal to noise ratio to obtain the ac- 
curacies of Fig. 1 is approximately 15 db and is met by 
100 watt (peak) beacon, which weighs 10 Ibs. and re- 
quires 25 watts of input power. 

The tracking radar measurements must be trans- 
formed from the tracking station coordinate system to 
the inertial coordinate system. The ground based digi- 
tal computer, which is used in conjunction with the 
tracking radar, can easily perform the matrix trans- 
formation. 

The orbital elements needed in the rendezvous proc- 
ess will be relayed to the computer in the “control 
center’ by a radio link. Transmission bandwidth 
could be conserved by storing the elements of the ex- 
pected orbit in the computer prior to launch, and send- 
ing corrections to these values after the actual orbit 
elements have been determined. 

The “‘control center’’ and the ‘‘rendezvous”’ satellites 
are identified by coded beacons in each. The FPS-16 
could determine the relative position of the two satel- 
lites to considerably less than one mile at the time of 


observation. The tracking accuracy of the FPS-16 is 





Position Accuracy : 
Azimuth Angle 0.1 mil r.m.s. 
Elevation Angle 0.1 mil r.m.s. 
Range 5.0 yards r.m.s 

Rate Accuracy With 20 Sec. Data Smoothing: 
Azimuth Rate 0.0025 mil/sec 
Elevation Rate 0.0025 mil/sec. 
Range Rate 0.3 ft./sec. 





Fic. 1. FPS-16 Tracking Accuracies. 





Information courtesy RCA 
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not sufficient, however, to effect the rendezvous maneu- 
ver itself. It is necessary then that some other instru- 
mentation be aboard the satellites to determine their 
relative position and velocity. A coordinate system 
must be established from which the position and ve- 
locity measurements can be referenced. 

It is assumed that the assembled satellite vehicle is 
always to be aligned to the earth’s local vertical. 
Hence, ‘the coordinate system is rotating and makes 
one revolution in one period of the satellite. 

The local vertical can be determined by infrared 
scanning the earth disc. The system is composed of 
four infrared telescopes equally spaced around the 
satellite body and in the plane perpendicular to the 
local vertical. Each telescope determines the edge of 
the earth disc by detecting the earth-space temperature 
gradient. The local vertical is established by bisecting 
the angle between each opposing pair of infrared tele- 
scopes. 

The satellite rotation will be nearly constant since 
the orbit is very close to being circular. Once the an- 
gular rotation is imparted to the satellite, only small 
angular rate changes will be required from the average 
value as detected by the horizon scanning system. 

The problem of locking onto the sun or moon can be 
eliminated since both the sun and moon have their 
peak radiation at wavelengths other than that of the 
earth. By proper filtering techniques, only the earth 
disc will be determined by the system. Since there is a 
particular angular relationship between the telescopes 
for tracking the earth from a given altitude, the local 
vertical information is ignored if the relationship is 
out of a specified tolerance. 

The infrared telescopes are limited in determining 
the earth disc to approximately 2 milliradians in the 
200-500 mile altitude region by anomalies of the atmos- 
phere and the surface. When the overall system is 
considered—-i.e., the telescopes, telescope mechanical 
alignment, the orientation control system——the local 
vertical can be determined to at least 0.5°. 

The control system to rotate the satellite about its 
center of mass is determined from the time the system 
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Fic. 2. Possible satellite assembly. 
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Fic. 8. Control center orientation system, typical 5,000-Ib 
satellite 


must be oriented, the degree of accuracy desired, the 
weight, and the electrical power available. 

The control systems for the ‘“‘control center’ and all 
of the succeeding ‘‘rendezvous’’ satellites utilize an 
infrared system to determine the earth’s local vertical, 
but different reaction control systems. The ‘“‘control 
center’ will be in orbit longer than any of the ‘‘rendez- 
vous” satellites and will have the responsibility of 
orienting the entire vehicle once assembled. The 
“control center’ control system utilizes reaction wheels 
to eliminate the large amount of reaction jet propulsion 
that would be required. The reaction wheel system 
would be augmented by reaction jets, but used only 
when the capacity of the reaction wheel system is ex- 
ceeded or to correct for external torques. 

The ‘‘rendezvous”’ satellites will require orientation 
only for approximately one day prior to its being joined 
with the ‘“‘control center.’’ Since the ‘“‘rendezvous’’ 
satellite’s payload is predominantly the scientific equip- 
ment for the missions requirements, sufficient electrical 
power will not be available for reaction wheels. For 
the short period of time of operation, the reaction jet 
better serves the requirements. 

A reaction jet system has been previously investi- 
gated which will provide orientation to the local vertical 
to an accuracy of approximately 0.5° which requires, 
for a typical 5,000-Ib. satellite, only 10 Ib. of propellant 
day. The necessary thrust is obtained from a hydrogen 
peroxide system. 

A particular axis of the ‘“‘control center’ will be 
chosen to be aligned to the velocity vector by the 
mounting axis of the roll rate gyro. Four “‘rendezvous’’ 
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satellites are to be attached to the “‘control center’’ for 
the mission of Fig. 2. 

The 1-3 axis of the control satellite is assumed to be 
determined from the roll rate gyro mounting axis. 
rhe output of the roll rate gyro will be zero only when 
the 1-3 axis is parallel to the velocity vector. The 
polarity of the roll rate signal is determined by observ- 
ing the polarity of the pitch rate signal. Although 
theoretically possible, the 180° yaw orientation is 
unstable and will be prevented by system noise. 
The “control center” is oriented to the proper yaw 
position (0°) to ‘‘take on’’ rendezvous satellite number 
l. 

The “‘rendezvous’”’ satellite is oriented to the velocity 
vector with the longitudinal axis along the velocity 
vector by its roll rate gyro error signal. 

After the first rendezvous satellite is attached to the 
“control center,’ the “control center’? must rotate in 
the yaw plane to orient position 2 along the velocity 
vector. The pitch and roll rate signals are interchanged 
and the alignment proceeds as previously discussed. 

To orient the “control satellite’ in yaw for “‘rendez- 
vous’ satellites numbers 3 and 4, the error signals of the 
roll and pitch rate gyros are reversed, respectively. 

If more than four satellites are to be coupled to the 
‘control satellite,’ the roll or pitch rate signals cannot 
accomplish the yaw orientation singularly. Any yaw 
orientation angle can be achieved by resolving the roll 
and pitch rate signals. By positioning the resolver to 
the desired yaw angle with respect to the velocity vec- 
tor, then the ‘“‘control center’’ will rotate in the yaw 
plane until the resolver output is nulled, at which time 
the vehicle will be at the desired yaw orientation. 

‘control center’’ is 


The orientation system for the 
shown in Fig. 3. This system utilizes reaction wheels, 
augmented by reaction jets. The ‘‘rendezvous”’ satel- 
lite orientation system is shown in Fig. 4. Due to the 
short time that this system is to be in operation, reaction 
jets are used. 
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Fic. 4. Rendezvous satellite orientation system, typical 5,000-Ib. 
satellite. 
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Fic. 5. Pulse position modulation and decoding 


The mission will dictate the attitude the “control 
center’ and ‘‘rendezvous”’ satellite are to assume. The 
“control center’ could be oriented with its longitudinal 
axis along the local vertical and the longitudinal axis of 
the ‘‘rendezvous”’ satellites oriented to the local hori- 
zontal so that the ‘‘control center’? would be the hub 
and the ‘‘rendezvous”’ satellites the spokes (Fig. 2). 

After each of the satellites are in orbit, they are 
stabilized in all three axes by nulling the output of the 
rate gyros (S-1 open, Figs. 3, 4). Once stabilized, the 
satellite is programed through a search sequence to 
acquire the earth. 

With S-1 closed, the earth disc will be sensed by the 
horizon scanning system and the satellite will be oriented 
in pitch and roll to the local vertical. The desired yaw 
angle is then effected by closing S—2. 


Satellite to Satellite Tracking 


Measurements of position and velocity are made from 
the ‘‘control center’ and referenced in its coordinate 
system. A small airborne type tracking radar is con- 
tained in the “control center’’ to measure the relative 
position of the ‘“‘rendezvous’’ satellite. The position 
information versus time is fed into the airborne com- 
puter to calculate the relative velocity. 

The radar, computer, and command generation equip- 
ment are contained in the ‘control center’ to effect 
rendezvous of all succeeding satellites. Therefore, the 
majority of the available payload of each ‘‘rendezvous 
satellite’ is available to satisfy the overall mission re- 
quirement. 

A beacon would be incorporated in each of the “‘ren- 
dezvous”’ satellites to identify them to the ‘control 
center’ so that they can be joined at the proper “‘re- 
ceptacle.’’ The beacon output power needed is less 
than 500 watts (peak). Beacons in this power range 
are readily available. 

The operating frequency of the airborne tracking 
radar and the beacon is approximately 10,000 Mc. This 
frequency choice is made since such radar presently exist 
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Fic. 6. Infrared homing guidance system 


and their antenna can be made quite small while provid- 
ing the desired narrow beam. To operate at ranges of 
200 miles with beacons presently available requires 
that the radar peak power be approximately 50 kw. 
The average input power requirement would be 150 
watts. 

The minimum range of the pulse radar is established 
by the pulse length. The airborne radar considered 
here has a pulse length of 0.25 usec. The minimum 
range derived from this pulse length alone is 106 ft. 
Since some time must be allowed for recovery time of 
the TR and ATR devices, a minimum range of 150 ft. 
has been set. 

From the orbit elements, determined by the ground- 
based computer, the launching time of the “‘rendezvous’’ 
satellite is defined. Due to timing, propulsion, and 
guidance errors, it is not possible to effect rendezvous 
directly. The ‘‘control center’’ is informed when the 
“rendezvous” satellite is to be launched and at what 
time the radar aboard should try to acquire the “‘ren- 
dezvous’”’ satellite. At the specified time, the radar ini- 
tiates a conical scan program to acquire the ‘‘rendez- 
vous’’ satellite. The radarwill lock on the’ rendezvous’”’ 
satellite prior to it being injected into orbit. Although 
the radar is locked on and the airborne computer is 
operating, the main stage propulsion is controlled by the 
“‘rendezvous”’ satellite’s booster guidance system. The 
radar measurements and the computer results control 
the injection and vernier propulsion of the rendezvous 
rocket. The command link is through pulse position 
modulation of the radar. 

The commands to the “‘rendezvous’’ satellite are the 
two angles which the thrust vector must be pointed 
from the local vertical of the ‘‘rendezvous’’ satellite 
and the duration (the thrust magnitude is fixed). 
Although the computations are made from the coor- 
dinate system of the ‘‘control center,’’ the thrust vector 
of the ‘“‘rendezvous’’ satellite will use the generated 
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angles from the ‘“‘rendezvous’’ coordinate system. 
This is possible since the two coordinate systems only 
vary by 5 milliradians in the worst case. As the two 
satellites come closer together, the difference in the co- 
ordinate system is reduced proportionately. 

The commands are relayed to the ‘‘rendezvous”’ satel- 
lite by the tracking radar aboard the ‘“‘control center.”’ 
Pulse position modulation of the radar relays the thrust 
angles and the desired impulse. Consider, for example, 
a four pulse code, as shown in Fig. 5. The difference 
between the first and second pulses determines the 
amount of time the fixed thrust is to be applied, the 
difference between the third and fourth pulses indicates 
the plane of the thrust angles, and the difference be- 
tween the second and fourth pulses specifies the magni- 
tude of the thrust angle. The time interval between 
pulses can be chosen such that an ambiguity will not 
exist. 

During the last few hundred (perhaps 200) feet of 
closing the radar will not be functioning and an infrared 
tracker on the ‘‘rendezvous”’ satellite will steer it to final 
closure. During this time it will be closing at about 10 
ft./sec. and will also get one programed retrothrust to 
slow it to 2 ft./sec. for the last few seconds. Pitch and 
yaw will be controlled by small verniers similar to (or 
perhaps the same as) the orientation stabilization sys- 
tem described above. The infrared system is shown in 


Fig. 6. 


Equations and Dynamics of the Rendezvous 


The assumption that the control satellite is in orbit 
in a circular orbit of radius 7) and moving with velocity 
Vo = VGM, ro is supplemented by the assumption that 
the rendezvous satellite can be shot into a transfer orbit 
having apogee at approximately 7» and arriving at 7 at 
a position near the control satellite, preferably ahead 
of it, but traveling at somewhat lower speed. The speed 
difference is to be made up by firing the injection pro- 
pulsion, thus ending with both satellites in the same 
orbit. It is our purpose to show how the two satellites 
can be brought together and to demonstrate the feasi- 
bility of all of the elements used. 
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Fic. 7. Coordinate system. 
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The process is as follows: When the rocket in the 
transfer orbit approaches to within a few hundred miles 
of the control satellite, that satellite’s radar will acquire 
the rocket and its computer will determine when injec- 
tion propulsion is to be used and what cutoff velocity 
should be used. The computations are made directly 
in the relative coordinate system and the thrust vector- 
ing of the rendezvous rocket is done also in that system. 
Cutoff of injection propulsion will occur as close as 
possible to a condition which will make collision occur 
at some future time. Immediately after shutdown of 
the injection propulsion a small or vernier (on-off) sys- 
tem is activated and a correction made, and such correc- 
tions are made several times before rendezvous. During 
this phase the radar on the control satellite tracks the 
rendezvous satellite and the control computer uses 
position versus time to generate velocities and also 
solves the equations of motion to see if intercept will 
occur as predicted. When the expected error gets suf- 
ficiently large a correction in velocity is made. As the 
rendezvous satellite approaches within a few thousand 
feet of the control, retrothrust sufficient to cut the 
velocity down to a closing speed of 10 ft./sec. is applied. 
Then as the range decreases to 300 ft. the computer 
calculates time to collision and the retrorockets are 
set to take out all but 2 ft./sec. at a few seconds (3 to 5) 
before impact. Inside of 300 ft. the radar loses track 
and an infrared beacon is displayed by the correct 
conical receptacle in the control satellite. The rendez- 
vous nose cone now activates (perhaps from 1,000 ft. or 
more out) an infrared seeker which senses correct 
direction and through a jet control system steers the 
nose cone into the conical receptacle, where it is me- 
chanically grasped and the umbilical fit is made. 

To effect this dynamical picture we must have the 
equations of motion of a body close to the control 
satellite and with not greatly different velocity and in 
coordinates originating from the control satellite. For 
this we have chosen a right-handed (x, y, z) system with 
y in the direction of the radius from earth’s center to the 
control and with z in the direction of the rotation vector 
for the control satellites—i.e., perpendicular to the plane 
of the control orbit and in the direction of its angular 
momentum. It turns out then that x lies along the 
tangent to the orbit and positive x is in the direction 
counter to the satellite’s direction of motion (see Fig. 7). 
In this coordinate system the equations of motion are 


9 


F,/m = & — Quy — u*x 
F,/m = 9 + 2ux — w*(y + 70) 
F,/m = 3 


F,, F,, and F, contain gravitational forces and perhaps 
thrust. The magnitude of the total gravitational force 


is 
—GMM ceartn)/ (x? + (y + 10)? + 2°] 


For values of x, y, 2 such as we are considering 
/ mad ~ - . 
(Vx? + y? + 2? < 200 mi.) we can expand the gravita- 


tional terms in powers of V x2 + y? + 27/ro and drop 
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PLANE OF CONTROL SATELLITE ORBIT—, 


PLANE OF RENDEZVOUS SATELLITE ORBIT——— 


Fic. 8. Rendezvous satellite plane inclination 


terms of order higher than (x? + y* + 3°)/ro?. When 
this is done both centripetal terms are cancelled but in 
the y equation a new term comes in which has a form 
similar to its centripetal term. The z equation also has 
one term of this sort. The equations now are 


T,/m = * — 2uy 
T,/m = ¥ + 2ux — 3w*y 
T,/m = 3 + ws 


Inspection shows that an integral of the motion is 
t 

f (7,4 + Typ + T.3)dt = 
t 


[(1/2m) (x? + 9? + 8°) + (1/2)(GMm/r*) (2? — 3y?)]}, 


f T,dt = [m(% — 2uy)|;, 
if 


0. For this case the 


Also 


Il 


Consider next 7, = 7, = T, 
solution is* 


2X, 1 of 
Ni) q 0 a 2X0 
y= — dy} cos ut + sin ut + | 4¥ — 

w w w 


2x 24 
4 “XG o : -/0 
x=2 — 3¥0} sin ut — cos wi + 


- ‘ 
) u 


a 
” Qn 
(6u yo — 3Xo)t + | Xo + — 


w 


n 


30 - 
= 2 cos ut + sin «t 


a“ 


where Xo, Yo, 20, Xo, Vo, 39 are initial values of position and 
velocity coordinates. 

The meaning of the z equation is simply that from the 
control satellite the inclined plane of the orbit of the 
rendezvous satellite makes it appear to perform simple 
harmonic motion in the z direction. (See Fig. 8.) 

dz So 


= 0 when tan ut 
dt UZo 


I| 


or Zmar = V 2" 
and the velocity at z = 0 is 
(2)mar = WV 2? + (44/w)? 


The z equation has a drift term which implies then 
that if % — 2wy ~ O then a drift superimposed on 


* Since delivery of this paper at the IAS June meeting in Los 
Angeles, it has come to our attention that these equations ap- 
pear in Space Technology, edited by H. S. Siefert, and published 
in the summer of 1959. 
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periodic x, y, and z motion will result. If the equations 
are arranged properly it can be seen that the motion in 
the x, y plane is of an ellipse with a moving center 
(translating uniformly) or a cyloid. Much can be 
learned directly from the equations. However, for our 
purposes here it is merely necessary to note that we can 
find initial velocities Xo, ¥, 29 such that for given xo, Yo, Zo 
and for a suitably chosen time to rendezvous, 7, inter- 
ception can be effected and, with proper retrothrust at 
the end, rendezvous. Ignoring the z motion for a while, 
we solve for initial velocities %, ¥) such that the ren- 
dezvous satellite passes through the origin at time r. 


These are: 


Xo Xo Sin Wr + yolGwr sin wr — 14(1 — cos wr) } 
Ww 3wr sin wr — 8(1 — cos wr) 
Vo 2x0(1 — cos wr) + yol4 sin wr — 3wr cos wr} 
w 3wr sin wr — S(1 — cos wr) 


[In the limit 7 — O these go %) —~ —(xo/r) and i —> 
—(yo/7r) as they should.| Having chosen 7 the control 
computer now calculates X%», ¥ and through a command 
link and the injection + vernier propulsion on the 
rendezvous rocket the correct conditions are established. 

As time passes the computer rechecks these equations, 
say at time 7’ (calculating %)(7’) and jo(7’) and by com- 
paring Xy(r’), %(7’) measured with the calculated ex- 
pected values, a correction command can be generated 
for the rendezvous rocket. Finally when the rocket 
nears the control satellite the measured velocities are 
used to compute the final retrothrust and the procedure 
described earlier is followed. 

The original choice of 7 will be made such that with 
some reasonably small correction in 3) (which is along 
with z) expected to be very small) z will be 0 in a time 
less than 7 so that the z motion may be stopped im- 
pulsively as z > 0. There will be then no further drift 
or oscillation in the z direction. Since the orbit planes 
should differ by less than 2° this is not expected to be 
any problem at all with respect to 7 and lacking further 
restriction 7 will be chosen so that the closing velocity 
on the average will be about 100 ft./sec. This will allow 
a retrothrust of about 500 Ibs. to bring the (5,000-Ib.) 
rocket to a 10-ft./sec. closing speed in approximately 30 
sec. (during which time it will travel about 1,500 ft.). 
For an initial separation of some 100 miles the closing 
process would then take about 90 min. or | revolution 
at the altitude we have chosen. In selecting the value of 
7 the computer makes an additional calculation of the 
total impulse needed to effect retrothrust and to allow 
for corrections along the path. If the impulse (pro- 
pellant) available to the rendezvous rocket were to be 
less than the computed amount needed, then the com- 
puter would reject this value of 7 and select another one, 
larger (corresponding say to 80 ft./sec. closure speed). 


Computer 


The digital computer to accomplish the above calcu- 
lations should weigh no more than 50 Ibs. It requires 
range and angle from the radar as a function of time. 
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COMPUTER 
Fic. 9. Computer block diagram. 


By coordinate conversion, differentiation, and smooth- 
ing it generates all of the necessary information to be 
used in calculations. One of the calculations it makes 
first is the sum of the momentum changes which must 
be effected if rendezvous closure is to proceed at ~100 
ft./sec. If this turns out too large then the computer 
uses SO ft./sec. and computes a new closing interval (7), 
and so on, until it gets a suitable 7. With this and using 
the equations above X%» and are calculated and com- 
mands sent to the rendezvous rocket. This is periodi- 
cally checked along the way and corrections made. 
When Vx? + y®? becomes approximately 2,000 ft. retro- 
thrust is applied bringing the closing speed down to 10 
ft./sec. at around 500 ft. Also a command is sent to the 
rendezvous satellite at 3 sec. before impact cutting the 
velocity to 2 ft./sec. Fig. 9 shows a rough sketch of 


this computer. 


Conclusion 


The first satellite, the ‘control center,’ is placed in 
the desired orbit. 

The ‘‘control cente1’’ is tracked with a system of in- 
strumentation radars, such as the RCA AN/FPS-16. 
The radar measurements are fed to a central digital 
computer where the elements of the ‘‘control center’s”’ 
orbit are determined. From the orbit elements, the 
ephemeris of the orbit is determined which defines the 
launching time of the second satellite, the ‘‘rendezvous”’ 
satellite. The elements of the “control center’ orbit 
are relayed to the computer aboard the ‘‘control 
center’’ for use in calculating the rendezvous maneuver. 

The “rendezvous” satellite is launched when the 
“control center’s’”’ orbit plane and the plane at ‘‘ren- 
dezvous”’ satellite injection will be coplanar. 

Due to timing, propulsion, and guidance errors, it is 
not possible to effect rendezvous directly. The ‘‘con- 
trol center’’ is told when the “‘rendezvous’’ satellite is 

(Continued on page 674) 
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A Method of Calculating Velocity 
Distribution for Turbulent Boundary Layers 
in Adverse Pressure Distributions 


EARL M. URAM* 
United Aircraft Corporation 


Summary 


A new method of calculating the behavior of turbulent bound 
ary layers in adverse pressure distributions is developed which 
permits direct determination of the velocity profile rather than 
the gross integral parameters normally used to infer the general 
character of the boundary layer. The method offers the sim- 
plicity of algebraic equations coupled with the use of charts 
rather than the laborious simultaneous solution of coupled differ- 
ential equations required by existing methods. The method also 
affords, for the first time, a means of determining the total bound- 
ary-layer thickness, thus allowing calculation of the absolute as 
well as the nondimensional velocity distribution. 

The velocity profile is considered to be composed of two re- 
gions—an inner region which is described by the Law of the Wall 
and an outer region which is described by a function depicting 
the deviation from that law. The deviation function involves 
two parameters which are uniquely dependent upon the skin- 
friction coefficient and a third parameter which, for practical 
purposes, can be considered a constant. Since the entire ve- 
locity distribution was found to be almost uniquely dependent 
upon the local skin friction, serious doubt is cast upon the gener 
ally accepted ‘“‘history concept’’ which considers the outer region 
of the boundary layer to be dependent on integrated upstream 
conditions. 

Agreement between experimental velocity distributions and 
those calculated by the method presented here is generally very 
The analysis and calculation procedures which are pre- 
two-dimensional, pseudo-two-dimen 


good. 
sented are applicable to 
sional,t and axisymmetric conical flows. 


Symbols 
a, b, c = coefficients of deviation function 
Cr = skin-friction coefficient, 2[U*/U,]? 
D = diameter 
H = boundary-layer shape factor, 6*/@ 
Zz = characteristic length, inches 
R, =: free-stream Reynolds number, U,x/v 
Roe = momentum thickness Reynolds number, l’)6/y 
U = temporal mean velocity at a point 
Uy? = skin-friction velocity, V1, p 
x = distance along stream direction or surface 
y = distance normal to stream direction or surface 
6 = total boundary-layer thickness 
6* = displacement thickness, ; {1 — (U/U,)]dy 


Received July 16, 1959. Revised and received December 2, 


1959. 

* Formerly, Research Engineer, Research Department. Now 
with Applied Mechanics Section, Research and Development, 
Electric Boat Division, General Dynamics Corporation. 

+ Herein defined as two-dimensional channel flows with rela- 
tively low aspect ratio for the wall on which the boundary layer 


is grown. 


n = nondimensional distance from surface, y/6 
é 
6 = momentum thickness, f [1 — (U/U,)|(U/U,)dy 
0 ‘ 
v = kinematic viscosity of fluid 
p = density of fluid 
Tw = wall shear stress 
o = deviation function 
Subscripts 
0 = initial value 
] = free-stream value 
Introduction 


OST ATTEMPTS to develop analytical methods for 
M predicting turbulent boundary-layer develop 
ment in adverse pressure gradients have involved the 
solution of coupled differential equations, which can 
give only gross information as to the character of the 
parameters, momentum 
A critical analysis of such 


flow, such as the integral 
thickness, and shape factor. 
methods! shows that the momentum thickness can be 
determined quite accurately but that no consistently 
reliable method exists for determining the shape param- 
eter. A more direct approach for describing turbulent 
boundary-layer development, of course, would be to 
determine the actual velocity distribution at each 
stage of growth. Knowledge of the velocity distribu- 
tion completely specifies all characteristics of the 
boundary layer required for practical applications 

e.g., determination of drag and pressure losses, heat 


transfer, etc. 
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Only very recently have methods other than the 
application of simple power laws been used to describe 
or predict the actual velocity distribution in the bound- 
ary layer. In references 2 and 3 the boundary layer 
is considered to be composed of an inner and outer re- 
gion, and it is shown® by the use of a normalizing pro- 
cedure that the outer region of turbulent boundary- 
layer velocity distributions form a near-universal family 
of curves. This same conclusion and a function de- 
scribing this family of curves were independently 
deduced in reference | from a detailed analysis of 
outer-region velocity distributions obtained from UAC 
Research Department tests and test data available in 
the literature. 

The object of the study reported here was to develop 
a calculation method based upon this two-region con- 
cept for determining the velocity distribution of a 
turbulent boundary layer under the influence of an 
adverse pressure distribution at any point of its 
development. 


Development of Boundary-Layer Calculation 
Procedure 
Boundary-Layer Velocity Distribution 
Fig. 1 is representative of a typical turbulent bound- 
ary-layer velocity distribution, which can be repre- 
sented as 
U : yU* . yU* 
— =F (> = o(n) + 5.6 log {= (1) 
c v v 
where 7 = y/6 and U* = Vt. p. The inner region 
is represented by the Law of the Wall, 
U/U* = 5.6 + 5.6 log (yU*/»); O<n < m (2) 
and the outer region is represented by 


(U, — U)/U* = (1) — om) - 

5.6 log 7; nm S<n<1 (3) 
where 7 is the point of initial deviation from Eq. (2). 
The function ¢(n) is defined as the deviation function 
and was found in reference 1 to be approximated by the 


expression 
$(n) = 5.6 + a(n — m0) + ber” — (A) 


where the coefficients a, b, and c are dependent upon 
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the conditions influencing the boundary-layer growth. 
For practical purposes, ¢(n) should be taken as equal 
to 5.6 for 7 less than mo. 

A study of the initial deviation points obtained from 
many experimental velocity profiles revealed that » 
decreases as a boundary layer grows in an adverse 
pressure gradient, and rarely exceeds the limits 0.08 < 
nm < 0.15. It was found, as shown in Fig. 2, that it is 
possible to correlate np) with the momentum thickness, 8. 
Fig. 3 illustrates the effect on a typical velocity profile 
of a variation in 7) between 0.08 and 0.13, and it is seen 
that serious errors will not be introduced into the cal- 
culation of a velocity profile in the event that a priori 
prediction of @ is not extremely accurate. The initial 
deviation point can therefore be obtained from 


no = 0.0851 — 0.05 log @ (5) 


for 6 expressed in inches. 
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Fic. 3. Effect of initial deviation point variation on velocity 
distribution. 
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Fic. 4. Deviation curves for adverse pressure distributions; 
pseudo-two-dimensional flow. 








$(n) 


FUNCTION, 


DEVIATION 


k 
ati 
dist 
The 
a, l 
use 
larg 
slig 
a sh 
sent 
vel 
ign¢ 
witl 
beh 
I 
to h 
A st 
coe 
be 1 
sho 
dev: 
frict 


that 
of t' 
was 
witl 
The 
thar 
gest 
disa 

a 
a at 
Eq. 
dist 


th. 
ual 


O1T1 


Cute? | N's 





CALCULATING VELOCITY DISTRIBUTION 661 





20-- 
[ REF. | 


CONFIGURATION WZ 













xX -IN SYMBOL 
8.3 ° 

i¢6r— 18.3 
32.5 
47.0 
59.0 
69.0 
79.0 


> (7) 
eoq@gsx &@ © GO 





FUNCTION, 


DEVIATION 
i 
je 





| | 
0 0.2 0.4 0.6 0.8 1.0 
NON - DIMENSIONAL DISTANCE FROM SURFACE, 7 


aud 





O it 


Fic. 5. Deviation curves for adverse pressure distributions; 


pseudo-two-dimensional flow. 


Figs. 4-7 demonstrate the typical agreement of devi- 
ation distributions obtained by the use of Eq. (4) with 
distributions for several experimental boundary layers. 
The correlation of Fig. 2 and values of the coefficients 
a, 6, and c obtained from the experimental data were 
used with Eq. (4). It may be seen that in cases of very 
large deviation rates the exponential term in Eq. (4) 
slightly overcompensates at about » = 0.9 and causes 
a slight dip in the curve in that region. This dip repre- 
sents usually less than 3 per cent error in the calculated 
velocity ratio at that point and, in practice, can be 
ignored by fairing the curve into the boundary value 
with zero slope or by using a higher value of c (the 
behavior of this coefficient is discussed in detail later). 

In order to calculate a velocity profile it is necessary 
to have means of determining the coefficients of Eq. (4). 
A study of the behavior of experimental values of these 
coefficients was made to find correlations which could 
be used in a calculation procedure. It was found, as 
shown in Fig. 8, that the boundary value ¢(1), of the 
deviation function can be correlated with c,, the skin- 
friction coefficient. It is seen from the expression 


(1) = 5.6 + a(1 — m) +3 (6) 


that the boundary value of ¢(n) is essentially the sum 
of two of the coefficients of the deviation function. It 
was also found that the coefficient a@ can be correlated 
with the skin-friction coefficient as shown in Fig. 9. 
The conical flow data‘ are seen to be consistently higher 
than the other data, so that a separate curve is sug- 
gested for use with such flows. The data of reference 8 
disagrees with all similar data in all correlations. 

The coefficient c is not as rigidly determinable as are 
a and b, in the sense that c must be calculated from 
Eq. (4) by using some chosen point in the experimental 
distribution and by using predetermined experimental 


values of the other coefficients. Study of the behavior 
of this coefficient revealed it to be sensitive to the point 
in the distribution chosen for its calculation. Syste- 
matic study of this coefficient revealed, as illustrated 
by Figs. 10-12, that there is no consistent variation in 
the behavior of the coefficient c for different types of 
boundary layers. The values of c shown in these 
Figures were computed from points in the ¢ distribu- 
tions lying between » = 0.74 and n = 0.85 selected 
so as to reduce the scatter in the c values. Fig. 13 
demonstrates the effect on a typical velocity distribu- 
tion of a variation in the coefficient c between the limits 
—0.15 and —0.30. In the light of these figures, it is 
reasonable to consider c as constant and equal to —0.20. 


Calculation Procedure 


In order to calculate the velocity distribution in a 
turbulent boundary layer under the influence of an 
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Fic. 6. Deviation curves for adverse pressure distributions; 
conical flow. 
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arbitrary adverse pressure distribution by using the 
information presented in the previous sections, it is 
first necessary to know or have some means of deter- 
mining the following: 

(1) Free-Stream Axial Pressure Distribution or Ve- 
locity Distribution—Depending upon the particular 
problem, this could be obtained either from experi- 
mental data or from ideal fluid flow considerations. If 
the latter is the case, it is possible to obtain more accu- 
rate results by an iterative procedure using the calcu- 
lated properties of the boundary layer for successive 
approximations. 

(2) Momentum Thickness or Boundary-Layer Velocity 
Distribution at Some Reference Station Within the Flow 
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Effect of variation of coefficient c on velocity 
distribution 


Fic. 13 


Field—In most instances, a calculation will start from 
the minimum-pressure point where it is reasonable to 
assume that the boundary layer closely resembles that 
developed on a flat plate at the same Reynolds number. 
A method of computing the properties of a flat-plate 
boundary layer is described in the Appendix. 

It is also necessary to have a method to obtain a 
reasonable approximation to the momentum thickness 
at the point of interest in the boundary-layer develop- 
ment so that m) can be determined from Eq. (5). Once 
the velocity profile is computed, it can be used to ob- 
tain a more accurate value of the momentum thickness. 

As pointed out and discussed in detail in reference 1, 
at least two methods of calculating the behavior of 
momentum thickness are available. For two-dimen- 
sional flows with strong adverse pressure gradients, the 
following relation? can be used: 


0/6 = [Uo U, |*-8 (7) 
For axisymmetric conical flows, Eq. (6) takes the form‘ 
@D ADo = [l ho l f |4.6 (7a) 


where D is the diameter of the conduit. 
Another method (cf. reference 1) of calculating two- 
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dimensional flows having no restrictions on the pressure 
distribution is to use the following equation: 


0.2155 x 
rr = | jon f Ut? dx + 
bah X¢ 


1/1, 2155 


(ons Us4), | (8) 


In reference 1 this relation was extended to apply to 
two-dimensional channel flows with relatively low aspect 
ratio for the wall on which the boundary layer is grown 
(herein termed pseudo-two-dimensional flows) and to 
axisymmetric conical and annular flows. The rela- 
tionship takes the form 
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jl p 0.2155 : x ne ; 
= VW U ieee 0.060 | l o 105 dx -+ 
4 * x 


\ 1/1.2155 


(6L.Rg-?!5 U5 105) } (9) 


where @ is expressed in inches and L is a characteristic 
length expressed in inches (such as the diameter for 
conical flows or the diffuser height for channel flows). 

(3) Local Skin-Friction Coefficient—In reference 10 
it is shown that, for two-dimensional and pseudo-two- 
dimensional flows, 


; 5.8 <. 10-5 if 6 ; 

“f= 0° .298 Ro , w < Cy * sata (10) 
: oA X& 1079 = 6 , 
me g!-41 BR,0.2785 ? isthe ie 1,000 (11) 


and, for conical flows, 


1.42 X 107% 

7 99.845 R,0.218 (12) 
These formulas yield skin-friction coefficients in agree- 
ment with those obtained from Law of the Wall calcu- 
lations applied to experimental velocity distributions. 
It is also shown in reference 10 that somewhat better 
agreement is possible for all flows with the modified Lud- 
weig-Tillman!* formula 


; 15.0 xX 10-* 8.05 \2.04 . 
er Rgo-1188 log 771.318 (13) 


It is possible to use Eq. (13) in the calculation pro- 
cedure. If the computation is started from a minimum- 
pressure point, an initial value of H corresponding to a 
constant-pressure boundary layer can be obtained as 
described in the Appendix. An alternate expression 
for calculating H for constant-pressure boundary layers 
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Fic. 17. Comparison of calculated and experimental velocity 
distributions; two-dimensional airfoil flow 


H = 1.876 — 0.134 log R,; dP/dx =0 (14) 


Since /7 is a slowly varying function of x, a step-by- 
step calculation of velocity profiles in the boundary 
layer can be made using the values of 7 from the initial 
velocity profile in Eq. (13) to determine the skin fric- 
tion of the next velocity profile to be calculated. If 
so desired, this process could be made iterative. 

Once the momentum-thickness approximation and 
skin-friction coefficient at a point of interest are deter- 
mined, the calculation of the velocity profile at that 
point proceeds as follows: 


(1) 70 is determined from Eq. (5) 

(2) (1) is determined from Fig. 8. 

(3) @ is determined from Fig. 9. 

(4) 6 is determined from Eq. (6). 

(5) cis assumed equal to —0.20. 

(6) The velocity profile is then calculated from 


U/Ur = 1 = [e,/2)[6(1) — 5.6 — a(n — m) — 
be~?°"-* _ 5 6 log n] (15) 


The correlation of ¢(1) with skin-friction coefficient 
(Fig. 8) can be used to obtain an equation describing 
the behavior of the boundary-layer thickness, 6. Eq. 
(1), expressed in terms of the boundary conditions at 
the outer edge of the boundary layer 7 = 1(y = 5), 
can be solved for 6 to yield the expression 


_ yp 2 0.5 - 1 J y 0.5 e i] 
aoe l= | 08" 56 = | o(1)¢ (16) 


Comparison of Calculated and Experimental 
Results 


A number of representative velocity distributions 
were calculated for the experiments of reference 1, 


4, 6, and 7, using the procedure described here. These 
distributions are compared with the experimental data 
in Figs. 14-19. Velocity profiles near the beginning, 
middle, and end of the pressure distributions of each 
of these experiments were chosen at random. Figs. 
14, 15, and 16 exhibit excellent agreement with the 
channel-wall data of references 1 and 6. The agree- 
ment with the airfoil data of reference 7 shown in Fig. 
17 is not as good, since those data account for some of 
the largest scatter in the correlations of Fig. 8 and in 
reference 10). 

The agreement shown in Fig. 18 with the conical 
flow data of reference 4 is considered acceptable, except 
for the distribution corresponding to x = 15.29 in., 
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Fic. 18. Comparison of calculated and experimental velocity 
distributions; conical flow. 
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Fic. 19. Comparison of calculated and experimental velocity 
distributions; conical flow. 
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Fic. 20. Comparison of calculated and experimental velocity 
distributions; two-dimensional airfoil flow. 


the difference being caused by the scatter in the corre- 
lation used in reference 10 to determine Eq. (12). Fig. 
19 shows the agreement obtained when that scatter 
is eliminated by using experimental values of skin- 
friction coefficient. Agreement comparable to that 
of Fig. 19 can probably be obtained by using the sug- 
gested procedure involving Eq. (13) for determining the 
skin-friction coefficient. 

It was found in reference 10 that Eqs. (10) and (11) 
for skin friction cannot be used for the reference 5 
data, so experimental values were used in the calcula- 
tion procedure; good agreement with experiment is 
exhibited in Fig. 20. Comparable agreement should 
be obtained if the suggested procedure for application 
of Eq. (13) to determine skin-friction coefficient is used. 

As indicated by the agreement with experiments, the 
accuracy of the proposed calculation method is gen- 
erally quite good. Improvement upon the proposed 
method of calculating the velocity profile could be 
made by developing (1) more accurate and universal 
means of determining the skin-friction coefficients, (2) 
an improved correlation for mo, (3) a method which 
accounts for the variations of the coefficient c, and (4) 
a means of correcting for the slight overcompensation 
of the exponential term at high rates of deviation from 
the Law of the Wall. The latter might be accomplished 
by enforcing the boundary condition d¢/dyn = 0 at 
7= 1. 


Appendix—Calculation of Turbulent Boundary 
Layers in Zero Pressure Gradient 


The zero pressure gradient turbulent boundary layer 
has been extensively studied in references 11, 12, and 
13. Such a boundary layer is self-preserving (cf. refer- 
ence | or 13); i.e., all of the velocity-defect profiles are 


SCIENCES—SEPTEMBER, 1960 


near-universal and can be approximated by a single 


curve. 

The local skin-friction coefficient at some position 
from the leading edge of a flat plate can be calculated 
from the Karman-Schoenherr equation (cf. reference | 
or 11): 


(cy)? = 1.67 + 4.17 log (R,c,), 10° < Rz < 109 
(17) 


Eq. (17) may also be used to approximate the skin fric- 
tion at the start of a pressure gradient if x is considered 
as the distance to the minimum-pressure point from 
the effective start of the turbulent boundary layer. 

An equation for calculating 6 can be obtzined by 
transforming the Karman-Schoenherr equation for 
flat-plate total drag coefficient to one for local skin- 
friction coefficient as a function of momentum thick 
ness Reynolds number. The result is 


Be 1 
= : (18) 
3 5.84 log [2U;6/v} 


A simpler but somewhat less accurate equation to use is 


v 0.0334 \-544 
log 6 = log —— + - 
U; * 


which is a rearrangement of the local skin-friction for- 
mula of reference 14. 

If cy and 6 are known, the velocity profile may be 
calculated by using the procedure outlined in the sec- 
tion entitled Development of Boundary Layer Calculation 
Procedure. However, it is possible to take advantage 
of the self-preserving nature of the flow and approxi- 
mate all zero pressure gradient profiles by a single pro- 
file. For flat-plate boundary layers, cy usually falls 
between 2.5 X 10-* and 5 X 10 Figs. 8 and 9 show 
that the parameters ¢(1) and a may be considered as 
constants in this range of c, with ¢(1) = 8.0 (ef. refer- 
ence 12) and a = 4.80. It may also be assumed that 
nm = 0.10 and is a constant. The velocity profile can 
then be calculated from 


U/U, = 1 — [e,/2]°5[2.4 — 4.80(n — 0.10) + 

1.82¢770"-* _ 5.6 log yn] (20) 
A more accurate momentum-thickness value can then be 
obtained from the calculated velocity profile. 
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An Extension of the Linearized Characteristics 
Method for Calculating the Supersonic Flow 
Around Elliptic Cones’ 


ANTHONY MARTELLUCCI* 
Polytechnic Institute of Brooklyn 


Summary 


The method of linearized characteristics as applied by Ferri to 
the flow about elliptic cones can be used to determine the surface 
pressure distribution, even when only linear terms are kept in the 
boundary conditions, provided an area rule requirement is satis- 
fied. In addition, the method can be applied for angles of attack 
provided the elliptic body geometry is specified in a manner that 
does not distort the cross section. The surface pressure dis- 
tribution obtained by this modified method is in reasonable agree- 
ment with experiment over the range of Mach numbers and 
semidiameter ratios considered. Experimental results for several 
conical bodies are presented 


Symbols 

a = speed of sound 

a/b = ratio of ellipse semidiameters 

Ap = base area of the cone 

# = pressure coefficient 

H = cone height 

= expansion factor of the ellipse 

VW = Mach number 

ry, ¥, 8 = spherical coordinates 

R = gas constant 

S = entropy 

v;, Up, W = velocity components in spherical coordinates re- 
ferred to the limiting velocity 

J = local resultant velocity referred to the limiting 
velocity 

a = angle of attack 

; = ratio of specific heats 

6 = inclination of the shock axis from free stream 

n = local inclination of the body surface defined by 
Eq. (10) 

y, = coefficients of the Fourier series defining the non- 
axially symmetric conical body surface 

y = coefficients of the Fourier series defining the non- 
axially symmetric shock surface 

Subscripts 

b = properties at the body surface 

n = properties of the superposed linearized flow fields 

Q = properties of the basic axially symmetric conical flow 

Ss = properties at the shock wave 
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= free-stream conditions 


= conditions on the conical body surface 


Introduction 


t p< TECHNIQUES FOR DETERMINING the flow past 
slender elliptic cones at zero angle of attack and 
low supersonic Mach numbers are rather well exploited. 
Interest has developed in predicting the flow properties 
for nonslender elliptic cones, inclined to the flow, at 
higher velocities. 

Most of the original theoretical investigations of the 
flow past slender elliptic bodies, such as the method of 
Kahane and Solarski,' proceed from Ward’s slender 
body theory which is based upon linearizing assumptions 
which are not valid at high Mach numbers. Even 
Van Dyke’s second-order slender body solution for an 
unyawed elliptic cone? cannot justifiably yield good 
results at high Mach numbers. 

The pressure distributions around conical bodies 
without axial symmetry can be obtained by the 
linearized characteristics method discussed in reference 
3. The application of this method to elliptic conical 
bodies has been considered in reference +. The present 
paper is an extension of the method presented in refer- 
ence + with special emphasis placed on extending the 
technique to large angles of attack. 

Basically, this method consists of the superposition 
of linearized flow solutions upon the solution of the 
nonlinear flow around a circular cone. The equations 
that define the velocity components of the linearized 
flow fields are solved by a step-by-step numerical 
procedure. This permits the determination of the 
shock wave created by non-axially symmetric bodies 
and, hence, gives solutions also at high Mach numbers. 
It was shown in reference 5 that satisfactory results can 
be obtained by the method of reference 4 for elliptic 
cones at zero angle of attack, even when only linear 
terms are retained in the boundary conditions, pro- 
vided that an area rule is satisfied. 

In the present paper, this method is extended to pre- 
dict the pressure distribution on conical bodies at 
angles of attack. This modified technique is applied 
to obtain the surface pressure distributions on several 
elliptic cones. Experimental results were obtained 
for the same cones at 1/4, = 3.09 and 6.0, and a com- 
with theory was Comparison was 


parison made. 
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Fic. 1. The polar coordinate system and velocity components, 


also made with several existing theories for determin- 
ing the pressure distribution on elliptic cones. 


Theoretical Analysis 


One technique for calculating the flow field around 
noncircular cones is the method of linearized character- 
istics.* This method proceeds, in general, by lineariz- 
ing the departure of the flow field of interest from some 
known nonuniform flow close to the actual one. The 
particular problem of the supersonic flow past a non- 
axially symmetric cone was treated by Ferri, by super- 
posing perturbations on the axisymmetric solution, the 
disturbance velocity components being expanded in 
Fourier series in the polar angle 3 up to the tenth term. 
The effect of small angle of attack within the accuracy 
of the linearized boundary conditions is considered by 
superposing upon the linearized flow solution of the 
non-axially symmetric cone at zero angle of attack, the 
flow solution for the equivalent circular cone at the 
angle of attack at which the actual body is being 
analyzed. With this approach the entire flow field 
and the shock wave can be determined over a wide 
range of Mach numbers, in contrast to orthodox meth- 
ods of linearization. For larger angles of attack, to rep- 
resent accurately the circular cone from which the 
elliptic cones are perturbated, the method outlined 
had to be modified, owing to the linearized boundary 
conditions. 

It was shown in reference 5 that satisfactory results 
can be obtained for the pressure distribution around 
elliptic cones for zero angle of attack from the method 
of linearized characteristics, even when only linear 
terms are retained in the boundary conditions, pro- 
vided that an area rule requirement is satisfied. In 
addition, the Fourier series expansion for the body and 
shock shape can be curtailed at the tenth term (i.e., 
cos 103) without any appreciable error. In reference 
5, a comparison was made for a specific body by re- 
taining terms up to cos 188; the resultant pressure 
coefficients are practically identical to those obtained 
by the approximation as far as cos 108. Thus the 
applicability of the method is based on the assumption 
of small perturbations and is related to the possibility 
of representing the body cross section correctly, rather 
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than to the problem of describing the perturbation flow 
field. 

The area-rule requirement for the expansion of the 
boundary conditions was also justified in reference 5. 
It requires that the zeroth term of the Fourier series 
representing the cross section of the body within the 
approximation of the linearized characteristics method 
should be obtained by an expansion satisfying the re- 
quirement of equivalent base areas. The Fourier 
body coefficients are defined as follows (Fig. 1): 


l - 
Wo, _ f v,(d)dd (1) 
T 0 
4 cs 
Vn = f vr(8) cos nddd (2) 
Tv 0 


where y(#) is defined from the geometry of the conical 
body to within an arbitrary parameter k. The value 
of k can be specified once the value of Wo, is selected so 
as to satisfy the condition 


m7 tan” Wo, = fA B (3) 


That is, the circular cone semi-apex angle yY,, from 
which the non-axially symmetric body is perturbed, is 
specified by the condition that the base areas of the 
elliptic and circular cones are to be identical. For the 
particular case of the elliptic cone at zero angle of 
attack, the parameter y,(#) defining the geometry of 
the ellipse is expressed as: 


f ka/H \ 
y,(3) = tan! - PET RESIS 
\[(a?/b?) sin? 8 + cos? 8}! 2f 

The parameter k can be obtained from the relation 
that results from substituting Eq. (4) into Eq. (1), 
once yY,, a, and 6 are selected. Filon’s numerical 
method of integration was used to solve the integral. 
Once the expansion factor k is known, the remaining 
body coefficients can then be determined from Eq. (2) 
explicitly. 

One possible representation of the flow field over an 
elliptic cone at angle of attack consists of the super- 
position of the yawed circular cone solution upon the 
unyawed elliptic cone solution. The resulting body and 





“az0® ae5* a=10° | 
a/b*1788 


Fic. 2. Body deformation resulting from the original Ferri 
method. 
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Fic. 3. Elliptic cone geometry for the angle of attack condition. 


shock shapes referred to wind axes are: 


10 


Yo -_ Wo, + p Vny cos no + @ COs v (5) 


n=2 


10 
¥; = Yo, + ) ¥ Y,, cos nv + 6 cos d (6) 
n=2 
At moderate angles of attack on the order of 5 to 
10 degrees, the above method does not yield suitable 
results. The basic cause of the discrepancy is the de- 
formation of the body. A calculation was made for 
an elliptic cone with a/b = 1.788 and H Vab = 3.315, 
which corresponds to the flow about a 16.783° semi- 
apex angle circular cone. The body shape that re- 
sults for this ellipse using Eq. (5) for 5° and 10° angles 
of attack is illustrated in Fig. 2. Thus the addition of 
the single term containing the angle of attack was not 
sufficient to express the ellipse at an angle of attack. 
Various attempts were made to “‘adjust’’ the elliptic 
body deformation that resulted due to the angle of 
attack. The method of approach was to retain the 
same form as Eq. (5) yet retain higher-order terms in a. 
Some forms that were used were: 


(Yoaxo = (Wo)a=-0 + acosd + 


a 
E N( Wn + Wnsi)y COS nd 
2Wo, x=1,3 
a? 
(Yoexo = (Wo)e-9 +t acosd? + cos 23 
" Op 
(Wo) ax0 = (Yo)a=0 + @ Cos ov a cos 29 + 
" Op 
a 11 
Aw cos nd + 
Wop n=1,3 
at 12 : 
a Q,, cos nd (7) 
4Yo>" n=0,2 
n— | n+ 1 
where \,, = = (Wr-1), — = n+1) 
(n — 2)? ] e 
Qn = 9 (Wn—2), = i 1) Wn + 


(n + 2)? — 1 
9 


- 


(Wn+2), 


In each case, the body that resulted was deformed from 
an elliptic shape. 

The problem was now reduced to one of geometry, 
that is, to develop a Fourier series representing an 
ellipse which is inclined at an angle relative to its cen- 


troidal axis (Fig. 3). This series is: 


10 
Y = Yo, + DY vn, cos nd (8) 
n=1,2 
where the coefficients y,,, are now functions of the angle 
of attack, and are defined by Eqs. (1) and (2), where 


, f(a b)* cos J sin a 
y,(0) = tan 1 


2 cos’ a 


(ka/H)? (a/b)* sin? & tan® a? |'/*) 9 
Q 2? ws 
a’ s 
and Q=1+ ( > — = ') cos? J 
b? cos? a 


The parameter k must be determined so as to satisfy 
Eq. (3) as before; however, now it is also a function of 
the angle of attack. In this manner, the body was 
found to remain elliptic to within the accuracy of the 
series selected for all angles of attack where y(#) is 
single valued and continuous. 

With the body defined in this manner, the method 
of linearized characteristics as presented by Ferri‘ 
was employed to develop the pressure distribution about 
elliptic cones for various angles of attack. 


Application of the Method 


An investigation was made using the present method; 
results were obtained at two Mach numbers, namely, 
M, = 3.08820 and 5.94090, for elliptic cones with an 
“equivalent radius’ H/Vab = 3.315. The ellipses 
considered had ratios of semidiameters a/b = 1.394, 
1.788, and 3.680. 
equivalent radius has a semivertex angle of 16.783 
Since tabulated flow field solutions were not available 
for this circular cone, both the nonlinear axisymmetric 
flow solution® and the linearized flow solution’ had to 
be determined for both Mach numbers. The coeffi- 
cients defining the cross sections of some of the various 
bodies considered and of the resulting shock shapes are 


The circular cone with the same 
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Pressure distribution on an elliptic cone with a/b = 
3.680 and M, = 5.94 
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Typical shadowgraphs of the flow for 4, = 3.09. 
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EXTENSION OF THE LINEARIZED 





a/b = 1.394, a = 15° 


tabulated in Table 1. It should be noted that the 
shock falls inside the body for the slender ellipse with 
a/b = 3.680. 

The pressure coefficient is obtained from the ex- 


pression 


2 1 — Ve2\777! AS 
Cp = a a ezp {| — — 1 (10) 
yM)° l— V;? R 


where | (v,)y,° + (wy 
g a ? 
» | _ ¢ 
and = (z,)y,7 = | (au, + do Wns(@r,)y, COS nd | 
i 
(w)y,” > ny,,(w)y, sin nd | 


n=1 + 
The entropy rise can be expressed in the form 


AS l ' 27M,’ sin? ¥. — (vy — 1) 
og, 


Roy] y+1 


| (y + 1)4/,2 sin? | 
y log, —e (11) 
(y — 1)M,? sin? y, + 2]! 


The pressure distributions determined by the present 
theory for the four cones are shown in Figs. + through 
7, for AN; 3.08820, and in Figs. 8, 9, and 10, for 
M, = 5.94090. 


Experimental Investigation 


Equipment and Technique 

The experimental investigation was conducted in the 
supersonic and hypersonic facilities of the Polytechnic 
Institute of Brooklyn Aerodynamics Laboratory. 

The supersonic facility’ consists of a high pressure 
air source supplying air to two blowdown wind tunnels. 
In the tunnel used for this investigation, a Mach num- 
ber of 3.09 was established by fixed, symmetrical, two- 
dimensional nozzle blocks, with a 10-in. by 10-in. test 
section. Tests were conducted at a stagnation pres- 
sure of 250 psia and a stagnation temperature of 500°R. 
The Reynolds number per inch, based on free-stream 
conditions in the test section, was 3.57 & 10° The 
running time per test was on the order of 25 sec. 
Several spark shadowgraphs, as well as photographs 
of the pressure recording system, were taken during 
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Typical shadowgraphs of the flow for 14, = 6.0. 


this time. A spark light source of approximately 
one microsecond duration was used for all tests. The 
hypersonic facility’ uses the same air supply as the 
supersonic facility. Air is directed from the air supply 
through a high temperature, high pressure heater" 
to an axially symmetric Mach 6 nozzle with a 12-in. 
diameter test section. The tests in the Mach 6 tunnel 
were conducted at a stagnation pressure of 600 psia 
and a stagnation temperature of approximately 1700°R. 


TABLE | 
Fourier Body and Shock Coefficients 
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The ratio of stagnation to wall temperatures was ap- 
proximately 2.3. Heat-transfer measurements for the 
elliptic cones have been made and will be reported on 
in the near future. The Reynolds number per inch, 
based on free-stream conditions in the test section, 
was 3.1 X 10°. The running time per test was on the 
order of 30 sec. Identical techniques and equipment 
for obtaining shadowgraphs and data were utilized. 

The models tested in both tunnels were elliptic in 
cross section. The ratios of the major to minor axes 
of the cross sections of the models tested were 1.00, 
1.394, 1.788, and 3.680. The ratio of the cone height 
H to the equivalent radius V ab for all the models was 
H/ Vab = 3.315 with the model height 7 = 4.54 in. 

The models were fabricated of 304 stainless steel of 
uniform thickness, and the static-pressure orifices 
were soldered to the model. These orifices were placed 
in one plane normal to the longitudinal axis of the body. 
Additional static orifices were placed along the } = 
180° and 3 = 270° rays to verify the conical nature of 
the flow. Stainless steel tubing of 0.062 in. diameter 
was used for the pressure leads. Depending on the 
pressure level to be measured, the surface static pres- 
sures were recorded by means of either 3-psi or 15-psi 
pressure differential transducers or a mercury manom- 
eter board. 
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Pressure distribution for an elliptic cone with a/b = 
1.788 and VW, = 3.09 


Fic. 14 


The model was mounted in both tunnels by means of 
a hollow sting attached to an angle sector. This angle 
sector permits a +10° range of angle of attack. In 
order to obtain flight attitudes above 10°, it was neces- 
sary to insert a 10° spacer at the aft end of the model. 


Results and Discussion 


The pressure distributions on elliptic cones with 
varying eccentricities (a/b = 1.00, 1.394, 1.788, 3.680) 
were determined experimentally at 1/; = 3.09 and 6.0 
at various angles of attack. Several typical shadow- 
graphs are shown in Figs. 11 and 12. The results ob- 
tained were compared to the pressure distribution cal- 
culated by the present theory. It can be seen from 
Figs. 4 through 9 that the theory agrees with experi- 
ment for each of the ellipses and Mach numbers con- 
sidered. For the slender ellipse (a/b = 3.680) the 
theory indicates that the shock falls inside the body, 
near the major axes. However, the pressure dis- 
tribution for the region where the solution is valid 
yields results that compare favorably with the experi- 
mental data. Hence, the ‘‘modified Ferri method’’ 
can be applied to cones whose cross sections differ from 
circular by as much as 3:1 with reasonable accuracy 
over a wide range of Mach numbers. 

Several theories exist that predict the pressure dis- 
tribution on elliptic cones for zero angle of attack. 
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The comparison is illustrated on Figs. 13 through 16. 
The method of Kahane and Solarski,! which proceeds 
from Ward’s slender-body theory, and Van Dyke's 
second-order theory,” which proceeds from slender- 
body theory and includes the effect of nonlinear terms 
in the equation of motion, cannot be expected to give 


good results at the speeds considered. 

The Newtonian theory for an elliptic conical body has 
been presented in reference 11. The pressure dis- 
tribution by this theory is given by 


C, 2 cos” 7 


where 7 is the local inclination of the surface— that is, 
(cos W)y~-9 sin asin y — (sin ¥)y-9 cos a V m* cos’ y + sin ¥ 19 
cos 7 = , = = a 
V m?K? cos? wy + sin® p 
where 
(sin ¥)5-0 (tan ¥)5=0 with experimental data for 1J = 6; however, at = 3 


i = m= 


(sin W) 5-90 (tan ¥)»=90 


The theory predicts pressures that are somewhat lower 
than experiment. 

The equivalent-cone method is presented in reference 
12. In this method an equivalent circular cone at 
zero angle of attack is defined for each 3 on the elliptic 
body, with the circular cone semi-apex angle defined 
by Eq. (10). 
the periphery of the elliptic body can then be obtained 


The surface pressure distribution around 


from reference 6 for each equivalent cone with the same 
free-stream velocity. The actual velocity components 
(v;, w) on the surface of the elliptic cone are not con- 
sidered. This method agrees surprisingly well for 
determining the pressure distribution on elliptic cones 
at \/, = 6.0; however, it is quite inconsistent at .1/; = 
3.09. 


Conclusions 


The modified method of linearized characteristics can 
be used to determine the flow about elliptic cones at 
angle of attack, provided the elliptic body geometry is 
specified in a manner that does not distort the cross 
section. Results obtained by this method agree with 
experimental results over the range of Mach numbers 
and semidiameter ratios considered. 

The results of the equivalent-cone technique agreed 
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the agreement was less consistent. The Newtonian 
approximation was consistently poor in agreement with 


experimental data. 
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Terminal Guidance System (Continued from page 658) 





launched and at what time the radar aboard should This procedure is repeated as necessary until the 
start trying to acquire the ‘‘rendezvous”’ satellite. The “rendezvous” satellite and the ‘‘control’’ satellite are 
only separated by 2,000 ft. At 2,000 ft. the vernier en- 


radar is used to acquire the “‘rendezvous”’ satellite prior 
gines are rotated 180° and retrothrust is applied so that 


to it being injected into orbit. The radar is locked onto 
satellite and the associated com- the satellites rendezvous and do not collide. The 
computer determines the amount of retrothrust which 
will reduce the closing speed to 10 ft./sec. To re- 


the ‘“‘rendezvous”’ 
puter is operating and the main injection propulsion is 
fired by the guidance system aboard the ‘‘control’”’ 


satellite. At the end of injection the vernier engines duce the initial velocity to 10 ft./sec., approximately 
are initiated and this propulsion unit controlled by the 1,500 ft. will be required. At the separation of 500 ft. 
radar aboard the ‘‘control center’’ makes the fine ad- a homing infrared system in the ‘‘rendezvous’”’ satellite 
justment in velocity. When considering the error of is signaled to acquire the ‘‘control center.” 
guidance and propulsion, plus allowing some freedom All commands to the ‘rendezvous satellite’ will be 
on the launch time, it appears that the two satellites from the command guidance systems until its minimum 
will be less than 200 miles separated. range of 150 ft. is reached. The closing is at 10 ft./sec. 
The satellite is attitude stabilized by an inertial as previously established, but it will be necessary to ap- 


ply side thrusts to rendezvous. When the minimum 


wheel/reaction jet system from information obtained 
range of the command guidance system is reached, the 


by an infrared horizontal scanner. 


The radar continually measures the relative position IR homing system takes over and the necessary side 
and velocity. From the initial measurement, the com- thrusts are then determined by this system. 
puter determined the thrust change that was necessary. At about 3 sec. before predicted impact time, the 
The following radar measurements are compared with closing speed is reduced to 2 ft./sec. by another burst 
the computed values for that time and, if different, the of retrothrust and final mechanical mating is ac- 
vernier (on-off) engines are fired to correct the error. complished by the remaining impact. 
++ + 
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A Sublaver Theory for Fluid Injection Into 
the Incompressible Turbulent Boundary 


Laver 


+ 


DONALD L. TURCOTTE* 


Cornell Unwersity 


Summary 


A sublayer region is introduced in which the intensity of turbu- 
lence grows at a prescribed rate. The decrease in wall shear 
stress due to fluid injection into the boundary layer is found under 
the hypothesis that the effect of injection is restricted to the 
Experimental measurements of the velocity 
The 


theoretical decrease in wall shear stress is in good agreement 


sublayer region 
profiles with fluid injection substantiate this hypothesis 


with experiment; the solution is particularly simple and for small 
values of the injection parameter it contains no arbitrary param 
eters. The theory provides a similarity parameter which differs 
from the one in general use. 


Symbols 


free dimensionless parameter 


skin-friction coefficient (27,/ ply? 


G = injection parameter 
u = velocity parallel to the wall 
Us = dimensionless velocity (u/u;) 
Ur friction velocity ( V 7,,/p) 
Nive friction velocity without injection (V 7,,,/p) 
U free-stream velocity 
velocity normal to the wall 
i injection velocity 
x coordinate parallel to the wall 
4 coordinate normal to the wall 
V4 = dimensionless distance from the wall (yu,;/v) 
€ = eddy viscosity 
p = density 
T shear stress 
Ta shear stress at the edge of the sublayer 
Tu shear stress at the wall 
T., = shear stress at the wall without injection 
v kinematic viscosity 


(1) Introduction 


ii IS WELL KNOWN that the injection of fluid into a 
flat-plate boundary layer will reduce the wall shear 
stress. For the laminar flow regime this reduction can 
be obtained by solving the Blasius boundary layer 
equation with a modified boundary condition to ac- 
count for the velocity at which fluid is injected. A 
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tabulation of the results has been given by Emmons 
and Leigh. ! 

Any solution to the analogous turbulent problem 
must include some degree of empiricism. Until there 
is a statistical theory of inhomogeneous turbulence, 
the solutions of turbulent boundary layer problems 
must rely upon the various phenomenological theories. 
The experimental logarithmic ‘“‘law of the wall’ ve- 
locity profile led Prandtl* to hypothesize a mixing 
length proportional to the distance from the wall. 
A solution for the decrease in wall shear stress due to 
fluid injection into a turbulent boundary layer has 
been given by Dorrance and Dore’ utilizing this mixing 
length hypothesis. Variations on this solution have 
been given by Clarke, Menkes, and Libby,‘ by Dor- 
rance,> and by Black and Sarnecki.® Since each 
solution contains at least one free parameter it is not 
surprising that some agreement with experiment is ob- 
tained. The principal objection to these theories is 
that they omit the effect of injection on the region near 
the wall, in which laminar exchange predominates. 

It was recognized by Prandtl that the mixing length 
hypothesis was not valid in the immediate vicinity of 
the wall. To account for this discrepancy he postu- 
lated a “‘laminar’’ sublayer in which turbulent exchange 
could be neglected. Concurrently with the work of 
Dorrance and Dore a theory including the effect of in- 
jection on a completely laminar sublayer was obtained 
by Rubesin.’ Again, reasonable agreement between 
theory and experimental shear measurements is noted; 
however, the arbitrary matching of the inner laminar 
solution to the outer turbulent solution provides several 
free parameters. 

It was recognized by von Karman* that in a number 
of problems even the two-region theory is inadequate. 
Therefore, von Karman introduced a third or buffer 
region in which both laminar and turbulent transfer 
are important. This somewhat 
modified by Rannie® who specifies an inner region in 
which the turbulent intensity grows at a prescribed 
rate, and an outer fully turbulent region. It is this 
two-region concept of the turbulent boundary layer 
which is applied to the injection problem in this paper. 
This choice leads to a shear law in which the value at 
the “outer edge’ of the sublayer region is reached 
asymptotically; as a result the solution is insensitive to 
the matching procedure. 


concept has been 
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(II) The Sublayer Region 


It is convenient to discuss turbulent boundary layers 
in terms of two parameters which are not actually de- 
fined physically. These are the mixing length, analo- 
gous to the mean free path in laminar momentum 
transfer, and the eddy viscosity, the turbulent analog 
of the kinematic viscosity. In Fig. 1 two alternative 
concepts of the flat-plate turbulent boundary layer are 
compared by plotting the dimensionless total viscosity 
(1 + ¢€/v) against the dimensionless distance from the 
wall (yu,/v). 

The value of the eddy viscosity corresponding to the 
Prandtl mixing length theory is 


e/v = 0.40, (1) 


The numerical constant is obtained by matching ex- 
perimental velocity profiles. In the Prandtl concept 
of the turbulent boundary layer, Eq. (1) is valid for 
y,> lland e/v = Ofor y, < 11. 

From experiments on heat transfer in pipes Rannie® 
obtained an empirical relation for the eddy viscosity 
very near the wall, 


e/v = sinh? (y,/13.89) (2) 


The numerical constant in this relation differs slightly 
from the value proposed by Rannie. This change is 
due to a different choice for the numerical constants 
in the logarithmic velocity profile. In the Rannie con- 
cept of the turbulent boundary layer, Eq. (1) is valid 
for y, > 25.4 and Eq. (2) is valid for y, < 25.4. 

Using the appropriate form of the momentum equa- 
tion, the velocity profile may be obtained once the eddy 
viscosity has been specified. The velocity profile 
corresponding to Eq. (1) is 
In y, + 5.10 (3) 


I 
u, = 
0.40 
Again the numerical constant of integration is obtained 
by matching experimental velocity profiles. The 
numerical values chosen are those given by Coles.!° 
The velocity profile corresponding to Eq. (2) is 
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Fic. 1. Dependence of the dimensionless apparent turbulent 
viscosity on the dimensionless distance from the wall without 
injection. 
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Fic. 2. Dependence of the dimensionless velocity on the dimen- 
sionless distance from the wall without injection 


These velocity profiles are compared with the experi- 
ments of Klebanoff'! and of Minkley and Davis’? in 
Fig. 2. 

It should be noted that a more sophisticated approach 
to the sublayer problem has been proposed. It can be 
shown theoretically '*'‘ that «/y ~ y,” for y* — 0, where 
niseither3or4. Relations for the eddy viscosity which 
fulfill this requirement aad which are asymptotic to 
Eq. (1) have been proposed by Reichardt'*® and van 
Driest.!° The complexity of these relations make their 
application to the present problem quite difficult. 
Since they are still empirical, the extent to which they 
are an improvement on the Rannie formulation is 
somewhat in doubt. 


(III) The Injection Solution 


The appropriate continuity and momentum equa- 
tions for a flat-plate incompressible turbulent boundary 
layer with fluid injection have been given by several 
authors.* +7 16 After integration these may be writ- 
ten 


vVv= Vw (5) 


T— puu = 7T, (6) 


And the shear stress is expressed in terms of the eddy 
viscosity as 


- ian 2. 
T= plv ae (7) 


It is hypothesized that Rannie’s relation for the eddy 
viscosity [Eq. (2)] remains valid with fluid injection. 
However, as the shear in the sublayer varies with in- 
jection, the appropriate value of the reference velocity 
u, is in doubt. For this reason Eq. (2) is modified 
to the form 


bys 


€ 
-= : h? = 8 
_ oe (+) (8) 


The factor 6 is a free parameter; its possible values 
will be discussed in detail when a solution is obtained. 
A substitution of Eqs. (5), (7) and (8) into Eq. (6) gives 
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SUBLAYER THEORY 
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sionless distance from the wall with injection 
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Introducing the dimensionless variables y, and uw, 
and the dimensionless injection parameter G = v,,/1,, 
Eq. (9) may be written 


he ( by, \ dit, G ; 
-osh? — Gu, = ( 
oi 13.89] dy, my 


The partial derivative has been replaced with the total, 
since no x dependence appears explicitly in the equa- 
tion. 

Equation (10) is a linear first-order differential 
equation which may be solved using an integrating 
factor, provided one boundary condition is specified. 
With the requirement that “, = 0 when y, = OQ, 
the solution is 


1 (13.89G/b) tanh (by, /13.89) 
Us = —le = 5 (11 
—: | | 


With fluid injection the shear stress in the sublayer 
region is, from Eqs. (7), (8) and (11), 
o = Tip 13.89G/b) tanh (by, /13.89) (12) 
For large values of y,, Eq. (12) has an asymptotic upper 
limit 
13.89G/b 
7. =t,0°"" (13) 


This asymptotic value will be identified with the shear 
stress at the ‘‘outer edge” of the sublayer region, which 
is a permissible conclusion if the asymptotic value 
(13) is approached within the region of validity of the 
original relation (8). According to the previous dis- 
cussion, the upper limit is y, = 25.4 without injection. 
If it is assumed that this upper limit does not decrease 
substantially with injection, then tanh (v,/13.89) 
is within 5 per cent of its asymptotic value. 

Now a second hypothesis will be made which was pre- 
viously proposed by Rott;'® namely, that the shear 
stress at the outer edge of the sublayer region is un- 
affected by fluid injection. The asymptotic shear 7,, 
and not Ty, is the same function of p, v, Up, and x with 
blowing as without. This assumption is equivalent 
to saying that the mechanism of turbulent momentum 
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transfer is so strong compared to laminar transfer that 
the effect of fluid injection will be limited to the region 
in which laminar transfer is important. This hypothe- 
sis will be substantiated when comparison is made with 
experiment. Since from Eq. (6) the value of the shear 
in the sublayer region without injection is a constant, 
the value of the shear stress at the outer edge of the sub- 
layer with fluid injection can be identified with the 
shear t,, which would have developed at the same 
station of the turbulent boundary layer without in- 
jection: 


Te = Tu (14) 


And from Egs. (13) and (14) the decrease in wall shear 
stress due to fluid injection is obtained as 
Tu/ Tu, _ es 89 (G/b) (15) 
Before a comparison can be made with experiment, 
the possible values of the parameter b must be deter- 
mined. With fluid injection the shear stress in the 
region falls in the range ry < 7 < Ty, From Eq. (8) the 
corresponding range in 6 is seen to be 1 < 6 < V tm Ties 
If the intensity of turbulence near the wall is unaffected 
by injection, then’ = V tm T», and the ratio of the wall 
shear stress with injection to the value without is 


— 13.89(vw/ uro) . 
Tw/Tm = eC ’ (16) 


If instead the turbulent intensity is determined by the 
value of the shear stress at the wall, then the factor } 
equals one. For this case, 


Te/ Te, = e —13.89(2w/uto) V two/ Tw (17) 


It may also be assumed that some ‘“‘true’’ intermediate 
value of 6 will be determined by some intermediate 
value of r. The assumption of the simplest possible 
average for b, namely 6 = (1/2)[1 + (u7o/ur) ], leads to 
the relation 


Te To, — e-°6 94(vw/uto) (1 + V rwo/Tu (18) 


However, for small values of the similarity parameter 
(v,/u,,) each of the above relations reduces to 
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Fic. 4. Dependence of the decrease in wall shear stress with 
injection on the dimensionless injection parameter. 
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fT v 
coi = pee— (19) 
T uy, 


Therefore this theory predicts the reduction in wall 
shear stress due to fluid injection without any arbi- 
trary parameters if the injection parameter is small. 


(IV) Comparison With Experiment 


Comparison of the above theory may be made di- 
rectly with the measurements of Mickley and Davis.'” 
These authors measured the velocity profiles of a tur- 
bulent boundary layer with injection through S0-mesh 
Jelliff Lectromesh screen. The corresponding values 
of the wall shear stress were obtained by using the 
appropriate von Karman integral momentum equa- 
tion. In Fig. 3 a set of the measured velocity profiles 
are plotted using the friction velocity based on the wall 
shear stress at the particular Reynolds number % ithout 
fluid injection. With these variables the profiles are 
clearly parallel. This substantiates the hypothesis 
that injection does not measurably affect the turbulent 
boundary layer outside the sublayer region. In the 
fully turbulent region of the boundary layer the Prandtl 
mixing length hypothesis is valid exactly as without 
injection. The displacement of the parallel profiles is 
due to the change in the velocity profile in the sub- 
layer, as discussed in Section (II). In principle, Eq. 
(11) could be used to determine the velocity profiles 
outside the sublayer; the arbitrary matching condition 
precludes this just as it does in the boundary layer 
without fluid injection. This property of the logarith- 
mic region with injection was first noted by Leadon 
and Bartle.'® 

The measured reduction in wall shear stress due to 
fluid injection is shown in Fig. 4. The experiments 
are compared with the present theory expressed in Eqs. 
(16), (17), and (18). The experimental points fall be- 
tween the limiting cases given by Eqs. (16) and (17), 
and Eq. (18) seems to be in excellent agreement with 


experiment. 


(V) The Similarity Parameter 


The present theory provides a similarity parameter, 
the ratio of the injection velocity to the friction ve- 
locity based on the equivalent shear stress without 
injection (v,,/u,,). This differs from the similarity 
parameter proposed by Rubisin” which is the ratio of 
the injection velocity multiplied by the free stream 
velocity to the square of the friction velocity (0,,0/,,”). 
Unfortunately it is difficult to distinguish between the 
two parameters on the basis of presently available 
experiments. A definitive experiment would have to 


discern the difference between c, and V cy. Because 
of the weak Reynolds number dependence such an 
experiment would have to be carried out over a Reyn- 
olds number range of 10° to 10°. The experiments 
of Minkley and Davis can be plotted using both simi- 
larity parameters. The results tend to favor the 
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present theory but are far from conclusive. 

However, a simple dimensional argument would 
seem to indicate the validity of the present theory 
over any parameter including a free-stream velocity. 
From the well-established applicability of the “law of 
the wall” for the region well out into the fully turbu- 
lent portion of the boundary layer, it seems clear 
that the only characteristic velocity in the turbulent 
injection problem is the friction velocity. Therefore 
the reduction in shear must be some function of the 
ratio of the injection to friction velocity, 


Tw/ Tu, = S(Vw/ UT) (20) 


This is just the result found from the present theory. 


(VI) Conclusions 


It seems clear that three aspects of the present inves- 
tigation indicate its validity: 

(1) The velocity profiles given in Fig. 3 indicate that 
the effect of injection on the turbulent boundary layer 
is restricted to the sublayer region. 

(2) The theory for small values of the injection 
parameter involves no arbitrary constants; it is in good 
agreement with experiment, as shown in Fig. 4. 

(3) For small values of the injection parameter only 
the region very near the wall should be affected; di- 
mensional arguments based on the characteristic ve- 
locity in this region (u,) substantiate the similarity 
parameter obtained from the theory. 

The validity of the present theory indicates that 
great care must be used to obtain experiments in which 
the sublayer region is not disturbed. This indicates 
a preference for measurements with an evaporating or 
sublimating material rather than the porous wall-type 
experiments available in the literature. 
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Analysis of Low-Aspect-Ratio 
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Summary 


Two methods are presented for the analysis of complex low 
aspect-ratio aircraft structures. Both methods provide for 
arbitrary external loading, are general with respect to the orien- 
tation of structural members, and permit arbitrary boundary 
conditions. For purposes of analysis a structure is idealized as a 
network of flexural members with interconnected torsion boxes 

In the first method, sets of linear equations are obtained by 
expressing boundary conditions, member deflection equations, 
equilibrium requirements, and slope-compatibility relationships 
in terms of deflections and internal forces. The solution for 
deflections and internal forces is then formed as the product of an 
inverse structural matrix and a column matrix of load functions 

In the second method, the conditions at a given boundary are 
assembled as a column matrix and are transferred in a step-by- 
step fashion over the entire structure to an opposite boundary. 
The transfer is accomplished by successive multiplications of 
square ratrices composed independently for the different trans 
fer ranges. The final operation is the inversion of a relatively 
small matrix and provides the solution for the unknown boundary 
conditions. 

Comparisons of theoretical results with experimental data and 


electric-analog solutions are favorable 


Introduction 


— TYPES of modern aircraft have low-aspect- 
ratio structural components. Such structures 
cannot be treated satisfactorily by the elementary beam 
theory; consequently, more refined analyses are re 
quired. The development of analytical procedures at 
the present, however, is restricted by two bounds. On 
the one hand, an analysis of deformations and stresses 
must have sufficient rigor to meet the demands of engi- 
neering accuracy; on the other hand, the capability of 
automatic computers places an upper limit upon the 
degree of complexity of usable procedures. 
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The literature of the past few years contains various 
methods for analyzing low-aspect-ratio aircraft struc- 
tures. It is not within the scope of this paper to pre- 
sent a detailed discussion of these methods; however, 
it will be well to note some of the more prominent con- 
tributions. In an early investigation, Schuerch! de- 
veloped a “wide-beam’”’ theory applicable to a series 
of beams and torque tubes. In the work of Stein and 
Sanders* at NACA, the wing is considered to be plate- 
like in behavior. Williams*® also developed a method 
in which he treats a wing structure as an equivalent 
plate. Levy® idealizes a wing structure as a network 
of crossing beams with interconnected torsion boxes and 
then develops a stiffness matrix for the complete struc- 
ture. Other investigators’—'’ also utilize the stiffness- 
matrix approach, but with differences in technique. 
Klein'': '* represents the entire structure by a system 
of force-carrying and shear-carrying elements. Equilib- 
rium equations and displacement-force compatibility 
relationships are expressed for each element; thus a 
matrix associated with both forces and displacements is 
produced. Using an idealized structure of flange and 
sheet elements, Argyris'* '* applies the unit-load and 
unit-displacement theorems to formulate a method of 
analysis completely in matrix language. In reference 
15 Argyris also provides a review of structural methods 
and an extensive list of references. 

The two present methods are directed toward pro- 
viding reasonable values for deformations and stresses 
of complex built-up low-aspect-ratio aircraft structures. 
Although the same idealized structure is used in both 
cases, the approaches are quite different. The first 
method is based upon the principle of consistent distor- 
tion; the entire assemblage of structural elements is 
considered and the inversion of a single structural ma- 
trix is required. The second method is derived from the 
concept of transfer matrices; segments of the structure 
are dealt with successively and the necessity of invert- 
ing high-order matrices is avoided. 


Structural Idealization 


The load-carrying capability of the structure is con- 
sidered to be provided by distinct flexural and torsional 
components. Spars and ribs are flexural members 
formed by assigning to webs, in addition to flanges 
which are directly attached, proportionate amounts of 
flange material from adjacent cover plates and stiffen- 
ers. Cover plates and webs also resist torsional shear 
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deformation. This characteristic is conveniently rep- 
resented by torsion boxes. Zero torsional stiffness is 
assumed for flexural members. A low-aspect-ratio air- 
craft structure suitable for analysis is illustrated in Fig. 
1. Solid lines indicate spars and ribs; dashed lines 
indicate torsion boxes which have their corners hinged 
to the intersections of flexural members. An enlarged 
view of a torsion box is given in Fig. 2. 

Before the stiffnesses of the idealized flexural and tor- 
sional components are established, the elements of the 
actual structure must be prorated. For example, the 
cross section of a typical built-up structure is shown in 
Fig. 3(a); the corresponding idealized cross section is 
given in Fig. 3(b). It is sometimes necessary to deter- 
mine the effective flange widths for flexural members 
having a relatively large spacing, or for a leading-edge 
member that crosses other members diagonally. Levy® 
suggests a treatment of these cases based on the effective 
width of wide beam flanges.'® This idealization as 
applied to a wing rib and a leading-edge member is 
shown in Figs. 4(a) and 4(b), respectively. In general, 
the effective flange width of a rib is a fraction of the rib 
length D_, but it should not extend on either side beyond 
the half-distance to the adjacent rib. For a leading- 
edge member, the effective width varies with the span 
length D, which is the span between intersections of 
spars with the leading-edge member. Full flange ef- 
fectivity is assumed for spars. The more haphazard 
the arrangement of spars and ribs, the greater is the 
need for intuitive judgment in idealizing the structure. 
Fortunately, structures encountered in practice usually 
have an orderly arrangement of spars and ribs. 

A representation of a fuselage structure as a network 
of beams and torsion boxes is given in Fig. 5. Because 
the fuselage itself is capable of resisting both bending 
and torsion, its action can be simulated by a series of 
torsion boxes and two flexural members located along 
the lines of intersection of the fuselage with the wing 
structure. 

In the mathematical model, spars and ribs are viewed 
as laterally unrestrained. Because of the plate-like 
nature of the structure as a whole, however, a certain 
amount of lateral restraint is exerted upon flanges of 
the flexural members. As pointed out in reference 6, 
if complete lateral restraint were provided, beam stiff- 
nesses should be increased by the ratio 1/(1 — v?), 
where v is Poisson’s ratio. Since only partial re- 
straints actually occur, a magnification factor of 
1/[1—(1/2)v?] seems to be a reasonable approximation. 

The torsion-box idealization permits the considera- 
tion of finite shear stiffness of all sides of the box. In 
an internal web, however, the shear flows from adjacent 
boxes tend to cancel one another and the best practice 
appears to be the assumption of rigid internal webs in 
the evaluation of torsion-box stiffness. External webs 
should be considered to have finite shear stiffness. 

Triangular torsion boxes sometimes occur along the 
leading edge of a wing structure. These boxes can 
usually be combined with neighboring quadrilateral 
boxes, as demonstrated in Fig. 1. 
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Consistent Distortion Method 


General Considerations 


In this method each straight-line flexural component 
establishes a member. For a symmetrical structure, 
the complexity of the analysis can be greatly reduced 
by consideration of symmetrical or antisymmetrical 
loading conditions, since any general loading can be ob- 
tained by a proper superposition of these separate con- 
ditions. It is convenient, then, to distinguish between 
two types of members. A ‘“‘through’’ member is one 
which has a line of symmetry coinciding with the line 
of symmetry of the structure as a whole; a “non- 
through’? member is any member which does not satisfy 
the requirement of a through member. 


Torsion-Box Reactions 
The corner forces applied to a general quadrilateral 
torsion box can be expressed in terms of the correspond- 


ing corner displacements as follows: 


7; Qi Qy Ay; Ais W 
T> - Q2, Gee Qe, A%4 2 1 ) 


T3 G3, G32 G33 Qs Ww 
Ta { G41 Ay a4 144 | W4 J 
where 
7T,...74 are the forces at corners 1, 2, 3, and 4, re- 


spectively, 
dy. ..@44, are terms dependent only upon the physi- 
cal properties of the box, 
and w;...%4 are the displacements of corners 1, 2, 3, 
and 4, respectively. 


An expression of this form, based on the work by Gar- 
vey” on quadrilateral shear panels, was derived by 
Thomann.’ Each corner of a box is considered to be 
attached to a particular member. It is possible, of 
course, that at a given point on a given member several 
boxes may contribute 7-forces. 


Member Deflection Equations 

If a flexural member is isolated as a free body, it is 
subjected to the following types of loadings: 

(1) Arbitrary external loading, which may consist 
of any combination of concentrated loads, distributed 
loads, and concentrated moments; 

(2) Cross-member reaction force, R;, applied to a 
member by a cross-member at point 7; 

(3) Torsion-box force, F;, the summation of all 7- 
forces applied to a given member at a point 7; 

(4) Joint moment, Q,, exerted on a given member by 
the joint at certain types of intersections. 

A typical member subjected to the various types of 
loads is illustrated in Fig. 6(a). Member coordinates 
are s and gz, as indicated. Two intersection points are 
selected as reference points for the member, and the 
origin of the coordinates is placed at the left-hand refer- 
ence point. The subscript / refers to the left-hand ref- 
erence point; the subscript 7 refers to the right-hand 
reference point. 
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Fic. 1. Typical structure for direct application of present 
methods. 





Fic. 2. Quadrilateral torsion box 
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(a) CROSS SECTION OF TYPICAL BUILT-UP WING STRUCTURE & ALLOCATION OF COVER SHEETS 


(b) CROSS SECTION OF (a) AS IDEALIZED FOR ANALYSIS 


Fic. 3. Typical built-up and idealized structure. 
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(a) RIB MEMBER b) LEADING-EDGE MEMBER 


Fic. 4. Effective flange widths of flexural members, suggested by 
Levy. 
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Fic. 5. Idealization of fuselage structure 
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Fic. 6. Member loads and displacements 
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Fic. 7. Deflected shapes of through member. 


' In responding to the various applied forces and 
i moments, a member is displaced to a deformed position 
as indicated in Fig. 6(b). It is convenient to consider 
the displacement to be a combination of a rigid-body 
translation and an elastic deformation. The two refer- 
ence points can be visualized as simple supports for the 
member free body, and the following geometrical rela- 
tionship can be stated for the rigid-body translation of 
a point 7 in terms of the translation of the two supports: 


Sn Sn 
(Wn) r = . a, oe l= a WW) (2) 


S, he 
The elastic deformation can be expressed 


(Wn) = DL Suki + Ze busFy + Do EuQe + mm (3) 


The quantities 6,,;, 6,;, and &,, are deflection influence 
coefficients and y, is an external-load function for the 
given member. These terms are defined as follows: 

(1) 6,; is the deflection at point m caused by a unit 
force at point 7; 

(2) 6,; is the deflection at point » caused by a unit 
force at point 7; 

(3) &,, is the deflection at point m caused by a unit 
moment at point k; 

(4) vy, is the deflection at point m caused by the ex- 
ternal loading. 

The final deflection at point 7 is now expressed as the 
sum of the right-hand sides of Eqs. (2) and (3). Thus 


Sy Sn 
Zz, = i WW, + (: ai *) U1 + be brik; + 
Yr “T t 
2 Or = x Ent Qe = Yn (4) 
J 
For a cantilevered member with the origin at the 


clamped end, Eq. (4) takes the form 
Ln = 2 bnilt; + » 5,jF; = »e EnkQ + Yn (5) 
a j 


If structural symmetry is used, Eq. (4) is modified for 


through members as follows: 


Sn Sn f 
wv, = WwW, a (: Pe uw, = pe 6ni*R; “F 
r S 1 


“rT 


De bnj*F;y + Lo Ene*Qe + Yn (6) 
J k 
Where the choice of signs exists, the positive sign applies 
to the symmetrical case and the negative sign applies 
to the antisymmetrical case. The asterisks have the 
following meaning: 

(1) For a symmetrical loading, influence coefficients 
are derived for two symmetrical unit forces (or mo- 
ments) ; 

(2) For an antisymmetrical loading, influence coef- 
ficients are derived for two antisymmetrical unit forces 
(or moments). It should be noted that in Eq. (6), is 
caused by the loading on the entire member. A through 
member distorted by a symmetrical loading and an 
antisymmetrical loading is shown in Figs. 7(a) and 
7(b). The deflection of the left support can be ex- 
pressed 


Ww; = +W, (7) 


where again the positive and negative signs apply to 
the symmetrical and antisymmetrical loadings, re- 
spectively. As Fig. 7 shows, corresponding R-forces, 
F-forces, and Q-moments at equidistant points on op- 
posite sides of the line of symmetry are equal in mag- 
nitude and have the same or opposite sense depending 
upon whether the loading is symmetrical or antisym- 
metrical. 

It can be seen from Eq. (6) that for a through mem- 
ber, deflections on half the member may be stated in 
terms of displacements and internal forces and moments 
located on the same half-member. If a member lies 
on the line of symmetry, one may use the artifice of 
“splitting” this single member into two members of 
half the actual stiffness spaced an infinitesimal distance 
apart. In this way consideration is still limited to half 
the structure. 

For every nonthrough member an equation of the 
form of Eq. (4) can be stated for each cross-point except 
reference supports. Similarly, for every through mem- 
ber an equation of the form of Eq. (6) can be expressed 
for each cross-point on the half-member being consid- 
ered, except the reference support. 

The influence coefficients and internal-loading terms 
are associated with statically determinate beams. 
Usually thin-webbed and large-flanged members are 
dealt with, and in these cases the incorporation of shear 
deformation is straightforward. The method of virtual 
work lends itself well to the calculation of these coef- 
ficients and is adaptable to numerical integration for 
digital computer operations. 


Member Equilibrium Equations 


Additional independent equations are obtained by 
expressing the equilibrium relationships for each mem- 
ber isolated as a free body. A summation of forces in 
the z direction produces 
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~ Ri + DF = -V (8) 


where V is the total applied external force in the posi- 
tive z direction. A summation of moments about the 
origin provides 
skit Ds F;+OGQ = -—M (9) 
; i 
where ./ is the counterclockwise moment of the total 
external loading about the origin. 

For through members, the condition of symmetry or 
antisymmetry was used in devising the special influence 
coefficients for pairs of unit forces or moments. From 
an examination of Figs. 7(a) and 7(b) it becomes appar- 
ent that only one independent equilibrium relationship 
exists for each through member. It is appropriate, 
then, to use the following equilibrium expressions: 

(1) For nonthrough members, both Egs. (8) and (9) 
provide independent equations; 

(2) For through members subjected to symmetrical 
loading, Eq. (8) should be applied considering only the 
half-member (zero shear occurs at s = s,/2); 

(3) For through members subjected to antisymmet- 
rical loading, Eq. (9) should be modified by replacing 
s; and s; by (s; — s,/2) and (s; — s,/2), respectively, 
and then applied considering only the half-member 


(zero moment occurs at s = s,/2). 


Boundary Conditions 

It is frequently convenient to use as a reference plane 
for a wing structure the plane passing through three 
points defined as follows: 

(1) A point of intersection of the longitudinal axis 
of the fuselage and a line parallel to the spanwise axis of 
the airplane; 

(2) The points of intersection of the trailing edge 
of the wing and the fuselage shell. 

Such a reference plane is established by the fictitious 
support system illustrated in Fig. 8(a). For conveni- 
ence, pin supports are used. For the accommodation 
of antisymmetrical loading, the single reference point 
By is connected to the fuselage shell by a rigid arm. 

A consideration of the boundary restraints in Figs. 
S(a), S(b), and S(c) leads to the following equations: 


up, = O for symmetrical and antisymmetri- 

cal loading (10) 
Up, = O for symmetrical loading (11) 
Rp, = 0 for antisymmetrical loading (12) 


For antisymmetrical loading, Rg, and Kz, must be zero 
in order to satisfy equilibrium requirements. 

It is, of course, possible to establish boundary condi- 
tions without resorting to the concept of pin supports 
and a rigid arm. It is helpful, however, to picture 
physical supports for the structure as long as such sup- 
ports do not alter the true nature of the structure. If 
the wing structure (with the immediately adjoining 
fuselage structure) of an airplane subjected to flight 


: y (b) RIGID ARM ~ SYMMETRICAL 
LOADING 
NGITUDINA . 


ANIS OF FUSELAGE = 


A 7 


(c) RIGID ARM ~ ANTISYMMETRICAL 
R R LOADING 
a’ SUPPORTS FOR ESTABLISHING REFERENCE PLANE 


Fic. 8. Plane of reference and fictitious supports 


loads is isolated as a free body, then appropriate shear 
forces, bending moments, and twisting moments must 
be applied at the “‘cut’’ sections where the forward and 
aft parts of the fuselage are separated from the wing 
structure. If these loads are applied as external loads 
to a structure with the fictitious supports of Fig. 8, the 


solution 
Re = Rp = Re = { 13) 


will be produced. In other words, these support points 
serve merely to establish a reference plane. When ex- 
ternal loads that are not self-equilibrating are applied 
to the wing structure alone, it is obvious that support 
reactions will be created. It might first appear that 
the effects of the external loads applied to the forward 
and aft parts of the fuselage are replaced by such reac- 
tions. One should note, however, that pin-support re- 
actions cannot represent the bending moments which 
ordinarily occur at the junctures. 


Joint Equilibrium and Slope-Compatibility Relationships 


The assumption of zero torsional stiffness for flexural 
members has been made. Therefore, if only two mem- 
bers intersect, neither exerts a bending moment on the 
other. However, when more than two members inter- 
sect in a rigid connection, the joint exerts a moment on 
each member. A member passing through such a 
joint, then, experiences a change in its bending moment. 
This induced moment is represented for a given member 
by the symbol Q,, where & identifies the point at which 
the change occurs. Additional unknowns are intro- 
duced by Q,-moments, and additional independent 
equations are needed. These can be obtained from 
equilibrium and slope-compatibility requirements at 
joints where such moments are introduced. If the 
angle y,, is the sweep angle of a member m, then at a 
joint k two equations of equilibrium can be stated: 


> Qem Sin Ym = 0 (14) 


m 


> Q, m cos Wm = 0 ( 15) 


m 


The slope ¢,m can be found in a manner similar to that 
used for deflections. Thus 





= 
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Fic. 9. Grouping of equations for final solution. 


nm = ie 
5 


$, S, 


w,+ > 0.2. + ¥ 6;F; + 
$ J 
: Bur Qk + Meg (16) 
k 


The quantities 6,;, 0,,, and 6,, are rotation influence 
coefficients for the given member; , is an external- 
load function. These terms are analogous to those as- 
sociated with member deflections and are defined as 
follows: 

(1) @,, is the rotation at point ” caused by a unit 
force at point 7; 

(2) @,, is the rotation at point ” caused bv a unit 
force at point j; 

(3) 8, is the rotation at point ” caused by a unit 
moment at point k; 

(4) X, is the rotation at point ” caused by the ex- 
ternal loading. 

These quantities can be determined by elementary 
beam theory. 

With the use of Eq. (7), Eq. (16) becomes for a 
through member 


l l 
Erm = ( ™ ) uy, “ig > Oni* Ry + » On3*F; + 
t J 


we s 


De Bue*Qx + Xn (17) 

k 
The asterisks have the same significance as those in Eq. 
(6). Pairs of unit loads or unit moments are applied 
symmetrically or antisymmetrically according to the 
symmetry or antisymmetry of the external loading. 
Where the choice of signs exists in Eq. (17), the nega- 
tive sign applies to the symmetrical case and the posi- 
tive sign applies to the antisymmetrical case. 

At a rigid joint the rotations of individual members 
may be related by geometric considerations. If ¢,; 
and ¢,, are the rotations at a joint in the x and y 
directions, respectively, the rotation of a member m 
having a sweep-angle yW,, is 


Ynm = ¥Ynz COS Ym = Pny sin Vm (18) 


A relationship in the form of Eq. (16) exists for each of 
N members entering the joint. Substitution of these 
expressions in Eq. (18) produces N equations in terms of 
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¢nz and ¢,,, as well as the quantities w, R, F, and Q. 
After the elimination of ¢,,; and ¢,,, (N — 2) equations 
remain. These equations, together with the joint equi 
librium equations (14) and (15), are sufficient to com 
pensate for the number of Q,-unknowns added to the 
problem. 


Grouping of Final Equations for Solution 


The equations collected for final solution are: 
(1) Member deflection equations [see Eqs. 
(4)—(6)] ; 

(2) Equilibrium equations [see Eqs. (8) and (9)]; 

(3) Boundary conditions [see Eqs. (10)—(12)]; 

(4) Joint equilibrium and slope-compatibility rela- 

tionships [see Eqs. (14)—(18)}. 

The diagram of Fig. 9 illustrates the grouping of the 
assembled equations. The basic unknowns appearing 
in these equations are w, RX, and Q. The F-forces are 
eliminated by substitution of their equivalents in terms 
of deflections. A collection of the external-load func- 
tions (y, V, M, and X) as right-hand sides of the equa- 
tions provides 


Tb by. ; a fw 7 ; & 7] 


boy bo». ai R = ( 19) 


: y 0 
L [eae —— 


The square matrix involves properties of the structure. 




















The column 
matrix of w-, R-, and Q-terms is the matrix of unknowns. 


It is designated the ‘“‘structural matrix.”’ 


The ['-matrix is composed of external-load functions. 
Premultiplication of both sides of Eq. (19) by the in- 
verse of the structural matrix yields 
yet. ty 
lo} 


where [db] —! is the inverse of the structural matrix. 


(20) 


It is convenient in preparing the various equations 
to identify separately forces and moments which occur 
at a single intersection on the different members. For 
the sake of simplification, however, one should observe 
that by the requirements of consistent distortion a single 
intersection will have a single deflection. Further- 
more, the interacting force between two members will 
have a single magnitude. 

From the solution of Eq. (20), the torsion-box corner 
forces are readily evaluated; the box analysis is then 
performed as a localized problem. All flexural members 
can also be treated as isolated structural components 
with known applied forces and moments (R,, /;, Q,, and 
external loading). 

Since deflection, force, and moment quantities ap- 
pear, the number of equations is relatively large. The 
matrix requiring inversion, however, contains many 
zero elements. Ordinarily several loading conditions 
are considered for a given structure. The matrix to 
be inverted is not dependent upon the nature of applied 
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loading, and once obtained, the inverse is available for 
use with any number of load matrices. An influence- 
coefficient matrix can be constructed for the complete 
structure from separate solutions for different unit 
loadings. 


Comments 

The method just described is similar in several re- 
spects to that of Levy.® Essentially the same struc- 
tural idealization is used, although in general torsion- 
box stiffnesses are evaluated somewhat differently. 
Also, in both procedures displacements and internal 
forces are related for each structural component and 
equilibrium relationships are satisfied. In the present 
method, induced joint moments are incorporated with 
member interaction forces and joint deflections to be- 
come the set of unknown quantities; the solution is 
then obtained by the inversion of a structural matrix, 
as previously discussed. In the Levy method, deflec- 
tions are the basic unknowns and are determined by 
the inversion of a stiffness matrix composed of smaller 
stiffness matrices for the various structural components. 


Transfer Matrix Method 


General Remarks 

The use of the transfer matrix concept avoids the 
construction and inversion of high-order matrices. 
The structure is partitioned into segments which are 
analyzed separately and independently. The re- 
assembly of the individual parts is accomplished by 
repeated matrix multiplications. This results in a 
matrix product relating the conditions at one boundary 
to those at another. The nature of this matrix equa- 
tion is such that the unknown boundary conditions can 
be determined by the solution of a relatively small set 
of simultaneous equations. In principle the method 
is not new and has been used by Thomson’? in the 
United States and Pestel,”? Marguerre,”?! Falk,?? and 
others in Germany. 


Application to Beam Structures{ 

The underlying idea in a structural analysis by trans- 
fer matrices is illustrated by considering the loaded 
beam of Fig. 10. At any section 7 of the beam the de- 
flection «;, the slope ¢;, the bending moment /;, and 
the shear force V; completely describe the condition, 
or state, of the beam at this section. The positive di- 
rections of these quantities are indicated in Fig. 11. 
Two cuts through the beam, one directly to the right of 
section 7 and another directly to the left of station 7 + 
1, isolate the 7th field of the beam. The four vector 
quantities w;,r, ¢:.2, Mir, and V;,z combine in a 
column matrix to form the state vector {Sir} at sta- 
tion 7,R: 

{ The application of the method to beams is presented here as a 


demonstration of the technique and is taken for the most part 
from references 20 and 22. 
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Fic. 10. Beam loaded by concentrated forces 
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Fic. 11. Coordinate system and sign convention 


- ’ a ¢ (21) 
(°a, RS M 
LV 


i,R 


Similarly, the state vector + Si4 LL} represents the quan- 
tities Wis1,2, Gi41,z, Migi,x, and Vi41,z. A linear rela- 
tionship between the state vectors }s;,z{ and }s,41,7} is 
readily established by the fundamentals of beam theory. 
In matrix notation this relationship is of the form (here 
neglecting shear deformation) 








a e 1,2 13 3 5 r _ 
7 2EI, GEL; ‘ 
en he 
a4 EI, 2EI; ’ 
V a ae 1! V 
it ha © 1c At ba 














where E is the modulus of elasticity and /; is the mo- 
ment of inertia in field 7; /; is the length of field 7. In 
more concise form, Eq. (22) may be stated 


BP i pie 2 - 
Sizt.z§ = [Pili Sirs (23) 


The square matrix [F;,], by which the state vector at 
station 7+1,L is related to that at station 7,R, is known 
as the ‘‘field”’ transfer matrix. The adding of the triv- 
ial equation 1 = 1 to the original set of four linear 
equations permits subsequent matrix operations. The 
length / of a field between two adjacent state vectors 
depends on the nature of the functions connecting them. 
External loading, intermediate supports, or certain 
internal conditions may require special matrices. The 
transfer of the state vector from station 1+1,L to 
station i+1,R, for example, is conveniently accom- 
plished by the introduction of a “‘point’”’ transfer matrix: 
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we | 7 Poy | ot a: ie , . 
kK: gE ! vi Saf = [Pal [Pra] [Fil... [Pel [lis 
¥ l Y 
Ml) = ! M (24) 9 OF 
= ; 
JE tied Ue {sa} = (LI{si} (29) 
= } oe L 1 Jett LI) Jette where [L] is the product of the field and point matrices, 
sited With this equation a direct matrix relationship is es- 
n short, . “4s 
tablished between the state vectors at the extremities 
fs. ' = ) or : se lalate ; 
Sittap = [Pisa] (S41.25 (25) of the structure. Explicitly stated, this connection 


If the bending stiffness of the beam of Fig. 10 is con- 
stant between successive stations, the relationship be- 
tween the state vectors at any two adjacent stations 
can be stated in the same manner as in Eqs. (22) and 
(24). Hence it is possible to obtain the following set 
of matrix equations: 


{ 52 u = [Fil tsi 

($2, ef a [P2] { s2,1} 

{ Siat,z} ee [Fil sir} (26) 
+ Sn} = [Fr—i)} Sn 1,R} 


These matrix equations are related in such a way that 
s j 

the internal state vectors {5»,,{ through }s,-1,,| 

can be eliminated. The substitution of 


io 2 gered. 2 97 
72,25 = [ Fi} ) say (27) 
in the second of Eqs. (26), for example, produces 


(sx,e} = [Po] [Fills (28) 


Further substitutions of intermediate state vectors lead 
to the equation 


{ { . 
between js,{ and {si} is 


Un = Ay + Apgi + A3My + ayVi + as 
En = AnWy + Agi + AM, + dogVi + os 
My, = zw, + a32¢1 + 33M, + d34V1 + G35) (30) 
Vn = Ag® + Ager + ay3M + a4gVi + C5 
l = l 


where the a,,-terms are known. The substitution « 
the boundary conditions 


w, = 0, fi = Ol 
w, = 0, o, = Of 


in Eq. (30) leads to two nonhomogeneous equations in 
MJ, and V;. The solutions for M4, and Vj, together 
with the known boundary conditions x; = 0 and ¢; = 0, 
establish the state vector at station 1. The state vector 
at station 2,L then becomes 


(. 1h) a 
{s2,cf = [Fill si} (32) 


Similarly, the state vector at any other station is ob~ 
tained from intermediate matrix products. 

It is seen that for the connection of the state vectors 
at stations 7,R and 7+1,R, the matrix product 
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is formed. In a more general sense the last column in 
Eq. (33) may be thought of as a load-coefficient column, 
since the effect of any type of intermediate loading 
would be reflected here. It is feasible, then, to transfer 
a state vector by means of a single matrix over a given 
range having an assortment of external loads. Also, 
the influence of a variation in the bending stiffness of a 
beam can be stated in terms of appropriate influence 
coefficients in the transfer matrix. This again would 
reduce the number of necessary matrix multiplications. 
In a continuous-beam problem, for example, transfer 
matrices for each span can be formulated directly with 
coefficients developed for the particular loading and 
bending stiffness of the span. 


Application to Rectangular Wings 


The transfer matrix method is also applicable to the 
analysis of multicell structures. The box beam of Fig. 
12, with idealized flexural members and torsion boxes 
of the sort already described, is given as an example. 

Since the torsion-box reactions are introduced as 
vertical shear forces at the intersections of spars and 


ribs, any field between two successive ribs may be 
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(b) VERTICAL FORCES ON SPAR-RIB INTERSECTION 


Fic. 13. Spar and joint free bodies 


visualized as a series of straight beam segments as 
shown in Fig. 13(a). The quantities w, ¢, WM, and V 
are introduced for each of the spars at the fuselage 
plane of symmetry and are transferred simultaneously 
over the entire span of the structure. The transfer 
matrices are of the order m = 4p + 1, where p denotes 
the number of spars. The transfer equations are the 
same for each spar and correspond to Eq. (22). In this 
example, however, the effect of beam shear deformation 
is included. The complete transfer matrix for a typi- 
cal field between ribs is given in Fig. 14. A trivial 
equation 1 = 1 is added to make the order of the field 
matrices conform to that of the point matrices. In an 
abbreviated notation this field matrix equation is 


> 


$141,L = F, SiR (34) 
The transfer of w, ¢, MW, and V from a station 
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Fic, 14. Typical field matrix. 
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i+1,L, immediately to the left of rib 7+1, to a station 
i+1,R, immediately to the right of rib 7+1, requires a 
new matrix. In rectangular wings with orthogonal 
spars and ribs the transfer equations 


nie = Witt | 
Giti.R = Pit1L ( (35) 
M41,R = M ist, 


apply if the displacements and the bending moments are 
continuous. The transfer of the shear forces across 
the rib, however, demands the incorporation of the ef- 
fects of (1) the external vertical loads, (2) the torsion- 
box reactions, and (3) the shear forces in the rib. 

The application of the equilibrium condition > F, = 0 
to a typical spar-rib intersection yields, according to the 
free-body diagram in Fig. 13(b), 


Viti V; — TP" 4+ T4547" - 
Tin + Ri + Pir (36) 


The 7-forces are reactions of the corresponding torsion 
boxes identified in Fig. 12. A,; denotes the rib force, 
and P;,; is the external load at spar JJ. Expressed in 
terms of rib influence coefficients and vertical displace- 
ments of the intersections, 


Ti+1,L 


Ry = [a2 A22 Ar 4] {uy Wy Wri u'ry} (37) 


Similar equations for the remaining spar-rib intersec- 
tions complete the set of transfer conditions at rib 7+1. 
Assembled in matrix form, the transfer conditions pro- 
duce the point matrix of Fig. 15. In short, this point- 
matrix equation is 


V; ay yz ays Ang 
Var az 422 423 au 
Vion a3) 432 a33 ase 
Viy as 242 das de4 
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Si+1 R Pig tI Jes S41 (38) 
;. +. t, 

0 l 
bigs t 1 


The torsion-box reactions are functions of the vertical 
displacements of the intersections of spars and ribs. 
For example, between spars J and // in Fig. 12 


7 I-II I-JT : , ‘ 
T; = k; (wr, i u7,) — (W754) — %'11;)| (39) 


— 7’ — 
i+2 W744) 


-— F-IT I-Il : 
Ti41 = Riss [(u', 
(U'717 115 — Wry | (40) 


It is evident that the reactions of the torsion boxes de- 
pend on the displacements at three successive stations 
1,7+1, and7+2. For the transfer of w, ¢, 7, and V 
across the rib, however, all equations must be stated in 
terms of transfer quantities at station 1+1,l. The 
unknown torsion-box reactions 

piu > 

fy) = | py (41) 
p fun 


are functions of the corner displacements at stations 

i+1 and 7: 

lr 

45 = 

[K |} u . . eh or. ¢@ . el (49 
LET, Bry Ey, SIV y4, B71; G77, 4711; GIV;5 \ 2) 


[AK] is the stiffness matrix of the three torsion boxes 
and is of the form 
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Fic. 15. Typical point matrix. 
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I-11 I-I1 
+k; _ k; 


[K,] = 4 piri _ pie 


I-IV II-1V; 
+ k; — k; 


The k-terms are the individual stiffness coefficients for 
the torsion boxes. 

Replacing the vertical displacements w in Eq. (42) 
by corresponding equations in Fig. 14 results in 


, Fi, 
ih = Rd [Ae | ted 


F;,, denotes the part of the field matrix in Fig. 14 that 
is associated with the deflection equations. The multi- 
plication of the rectangular matrices in Eq. (44) yields 
a new matrix [G]. Thus 


it, [G| + Semi (45) 


Substitution of the known relationships V;,zg = 
Viot,t, Mie = Misst — Vegi cli and 9,2 = I (Gisi.L, 
Misi.z, Visi,t) in Eq. (45) provides 


lef = [G*]\ Si41,2} (46) 
where [G*] is the appropriate modification of [G]. The 
matrix multiplication 

[Ji] tis _ [Ji] [G*] {Si LL} 


now permits the reduction of the matrix equation (38) 


Si+1,R (P; + J G*)'\ Figs Si+1,L 
— = cevaeeccefeeneme | | onsen (47) 


bras 0 | teat 
Before the torsion-box reactions 


tss3} nail [K is1] x 


os 9 op) 90) of" 
UF Tey0 S49 OT 49 SIVi42 144 


Cry Wry WV ay) 

(48) 
can be eliminated, a relationship must be established 
between the vertical displacements at stations 1+2 
andi+1. This is accomplished by premultiplying the 
point matrix of Eq. (47) by the deflection equations of 
the following field matrix [/;,;]. Three zero-columns 
are added to the original field matrix for the sake of 
conformability. The substitution of the resulting 
equations in Eq. (48) produces 


> Fizis(P; + JG*) a ree 
(teas = Kul] oe ee (Si4a,s tesa} 


(49) 
After simplification, Eq. (49) becomes 
5 tesat - [A ‘B)} Sisa,z, tisst 


where [| A :5] represents the product of the rectangular 
matrices. Thus 


tesa} = [A] }Si41.zf + [B] i tess} 


or (7) — (Bl) {tea} = [A] tsezf (50) 


i—pi 4 pri 
— pit 4 peu (43) 
— pv 4 pairs 





where [/] is a unit matrix. From Eq. (50) 
t 


ttesas. = (U] — (B))— TA] fsenacf = [D] bseer.c} 


(51) 


where [D] is the result of the indicated matrix opera- 
tions. The torsion-box reactions may now be expressed 
as functions of the known transfer quantities alone. 
By Eg. (51), 


[J i+] itis = [J s+] [D] Seat Lt (52) 


Substitution of Eq. (52) in Eq. (47) gives 


Si4i,R = (Pisa + JiG* + JigsD) Sit, x] 


(53) 

After the point matrices are reduced to the form 
shown in Eq. (53), both the field and point matrices 
are of the order = 4p + 1. The connection between 
the state vector at the fuselage plane of symmetry and 
the state vector at the wing tip is now established with 


the matrix product 


Serj = [Pal [Fr-al [Pr-il [Po] [Fil tsi.es 


or Snr} = [L] (ss R} (D4) 


From the insertion of the known boundary conditions 
M = Oand V = Oin }s,,2} and w = Oand ¢ = Oin 
+ Si.Rt> a set of simultaneous equations of order 2p 
results. The solution of this system determines the 
unknown terms in the state vector }s),z{. The multi- 
plication of the individual matrix products by this 
state vector finally yields the values of w, ¢, M, and V 
for all spars at each station. 


Application to Nonrectangular Wings 

The analysis of a swept- or delta-wing structure can 
be performed in much the same way as for the rectangu- 
lar wing already discussed. The procedure permits 
various choices with respect to the direction of the 
transfer process, the establishment of boundary and 
continuity conditions, and the construction of the in- 
dividual matrices. 

The treatment of certain problems which frequently 
arise with nonrectangular wings may be illustrated by 
a consideration of the wing of Fig. 16. The quantities 
w, ¢, M, and V for each spar are again introduced at 
the fuselage plane of symmetry and transferred in the 
direction of the wing tip. The leading edge is cut at 
the locations indicated in Fig. 16. The two straight 
spars thus created are treated in the usual manner. 
For the three kinked spars, the transfer process is car- 
ried across the kinks and then terminated at the cuts. 
Of the four quantities w, ¢, WM, and V which exist at 
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each of the cut faces, the bending moment / and the 
shear force V are incorporated into the transfer condi- 
In Fig. 17 are shown the 
The transfer 


tions for the adjoining spar. 
required equations at a typical spar kink. 
f the displacement w and the slope ¢ is discontinued; 
thus two unknowns are left at each cut. For their 
evaluation the conditions of continuity for displace- 
ments and slopes along the leading edge yield two 
independent equations at each kink. Beyond the 
wing tip the boundary conditions of zero slope and 
zero shear force lead to four more equations. Com- 
bined, these ten linear equations equal in number the 
ten unknown boundary conditions at the fuselage plane 
of symmetry. The unknown boundary conditions at 
the fuselage plane of symmetry are determined from 
these equations. The complete sets of boundary and 
continuity conditions for symmetrical and antisymmet- 
rical external loadings are given in Fig. 16. 

In a manner similar to the treatment of continuous 
beams, it appears feasible to analyze major aircraft 
components separately and then to combine the re- 
sulting matrix equations by a transfer procedure. The 
analysis of an entire aircraft could thus be accomplished 
with the use of relatively small matrices. 


Application of Methods* 


Deflections for the wing structure of Fig. 18 were de- 
termined by Benscoter and MacNeal?’ with the use of 
an electric-analog computer. The results of an analysis 
of the same structure by the first method (consistent 
distortion method) are illustrated in Figs. 19 and 20. 
The deflections of the tip rib for a concentrated load 
applied at different points along the tip rib are given in 
Fig. 19; the deflections caused by a pure twist loading 
are shown in Fig. 20. All loads are symmetrical with 
respect to the center line of the structure. 

A plastic-model test of the wing of Fig. 18 was re- 
ported by Zender** at the NACA. The available re- 
sults from this test are given in Figs. 19 and 20. 


* The consistent-distortion solutions given in this section were 
btained under the direction of Mr. F. J. Stebbins, of Convair, 
Fort Worth; the general procedure was programmed for the IBM 
704 digital computer by Mr. W. P. Durbin, also of Convair, 
Fort Worth. Appreciation is expressed for this help 


euemaw [ene eTHON 
sveenscens HEORY (SECOND METHOD - 
‘ 


ANALOG (REF. 25 





SPAR 
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Fig. 21. 
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Fic. 23. Delta wing with rectangular cross section (reference 
26). 


The electric-analog technique was also applied by 
Benscoter and MacNeal to the delta-wing structure of 
Fig. 21. A comparison of rib and leading-edge deflec- 
tions obtained from the analog and the second method 
(transfer matrix method) is given in Fig. 22 for the 
concentrated load applied symmetrically to each tip 
rib. 

A built-up metal wing with a delta planform was 
constructed and tested at the NACA.** The actual 
structure is shown in Fig. 23(a); the structure as 
lumped for analysis by the first and second methods is 
given in Figs. 23(b) and 23(c), respectively. Deflec- 
tions of the leading edge and ribs 9 and 11 as determined 
by test and by both theoretical methods are shown in 
Fig. 24 for a concentrated load applied symmetrically to 
each wing tip. 

As may be seen from Fig. 23, ribs are lumped in the 
structural idealizations of Figs. 23(b) and 23(c). The 
problem of accounting for torsional stiffness in the tip 
region arises. For the structure of Fig. 23(b), this 
stiffness is incorporated into a single quadrilateral box 
surrounded by the idealized flexural members at the 
tip. For the structure of Fig. 23(c), the tip region is 
not included in a torsion box; however, the flexural 
stiffnesses of the two tip members outboard of the last 
lumped rib are considered to remain constant rather 
than to diminish to extreme-tip values. 

A structural layout of the B-58 3/8-scale model is 








Deflection of ribs and leading-edge member for tip load 


Fic. 24. 
(symmetrical); wing of Fig. 23. 
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~ TABLE 1 
Influence Coefficients (Modified) for Deflection for Wing of Fig. 25—Analysis by First Method 


Deflection _—_ 
point 1 2 
1 1.00000 0 
2 0.85174 0 
3 0.72134 
4 0.53966 
5 0.28699 
6 





0.07109 


given in Fig. 25. The dashed line indicates a continua- 
tion of the line of rib 1 along the chord. This model is 
an example of a complex built-up wing structure. 
Experimental deflections for a concentrated load applied 
at point 3 (symmetrical loading) are given in Fig. 26 for 
various points along ribs 1 and 2. Analytical deflec- 
tions obtained by the first method are also given. The 
vertical coordinate is the ratio of the deflection at a point 
to the test deflection at point 3. The horizontal co- 
ordinate is the ratio C/C, (distance to deflected point/ 
root length); C and C, are identified in Fig. 25. 

Analytical solutions for concentrated loads applied 
at points 1, 2, 4, 5, and 6 were obtained in addition to 
that for a load at point 3. The modified influence- 
coefficient matrix of Table 1 gives the results of the six 
analyses. All terms result from dividing the actual 
influence coefficients by the deflection at point 1 caused 
by a unit load at point 1. The symmetry of the matrix 
gives an indication of the numerical precision of the 
computer operations. The matrix of unknowns for 
this problem contains 529 quantities. These include 
246 deflections, 229 internal forces, and 54 internal 
moments. The inversion of the structural matrix with 
single-precision machine operations required about 15 
hours. 

The curve of Fig. 27 provides a plot of order of ma- 
trix versus time of single-precision inversion of first- 
method structural matrices on the IBM 704 digital 
computer. The data are based upon analyses of struc- 
tures of varying degrees of complexity. The curve is 
thickened to show time-reductions obtained through 
refinement of the computer program. It is noteworthy 
that while approximately one hour is required for the 
inversion of a 240-order matrix, about ten hours are 
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Fic. 25. B-58 %/s-scale model (reference 27). 





——___—___-_-_-—— Loading point——____— 





0. 
0. 
0. 
0. 





4 5 6 
72136 0.53973 0.28700 0.07113 
65575 0.55842 0.26864 0.08412 
59621 0.44625 0.25687 0.06562 
44620 0.42618 0.20396 0.08531 
25686 0.20398 0.17563 0.04587 
06560 0.08528 0.04587 0.09627 


needed for the inversion of a 480-order matrix. Print- 
out times are not included in the times shown. 

The effect of shear deformation is included in the 
analytical results given. However, the consistent dis- 
tortion analysis of the structures of Figs. 23 and 25 
excludes the effect of flange taper on shear deformation. 


Conclusions 


The two present methods, while involving the same 
type of structural idealization, offer versatility in at- 
tacking problems of deflection and stress analysis of 
low-aspect-ratio aircraft structures. The consistent 
distortion method provides an approach which is rela- 
tively simple in application and which provides a ‘‘feel”’ 
for the structural behavior. Up to the limit of the 
capability of automatic computers to invert large-order 
matrices, this procedure permits a direct solution for 
deflections and internal loads. The transfer matrix 
method offers a means of avoiding the problem of large 
matrix inversions. It also provides a convenient way of 
combining analytically the behavior of smaller sub- 
structures within a larger structure. 

Other attractive features of the two methods include 
the following: 

(1) Both procedures lend themselves well to pro- 
gramming for digital computers. This is especially 
true for the transfer matrix method, since ordinarily 
many repetitive matrix operations occur. 

(2) An arbitrary distribution of applied loading is 
readily accommodated. 

(3) An arbitrary variation of stiffness along a flex- 
ural member is provided for. 

(4) No unrealistic restrictions are imposed upon 
root support or deformation of spars and ribs. 








07 08 


Cr 
Fic. 26. Deflection of ribs by a symmetrical tip loading for struc- 
ture of Fig. 25. 
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Fic. 27. Plot of order versus time for single-precision matrix inversion 


(5) Deformations and internal loads are considered 
simultaneously. Consequently, all structural com- 
ponents can be isolated as free bodies subjected to loads 
obtained directly from the primary solution. In other 
words, the determination of internal stresses is ac- 
complished by a consideration of applied forces and 
moments and not by a differentiation of deflection ex- 


pressions. 
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Theoretical and Experimental Investigation 


of Second-Order Supersonic Wing-Body 


T 


Interference 
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Summary 


Approximate second-order solutions for the supersonic flow 
around wing-body combinations are calculated, using two differ- 
ent theoretical models. In the first, the spanwise curvature of 
the body field is assumed small and the wing sweep small in com- 
parison with that of the Mach cone. In the second, two per- 
pendicular intersecting two-dimensional fields are considered. 
The analysis is restricted to such high Mach numbers that 17-?*< 
1, and an approximate formula common to the two models is 
then found for the second-order interference term. This for- 
mula can also be used to correct experimental pressure distribu- 
tions for the effect of nonuniformities in the wind-tunnel flow. 

In order to test the theory, wind-tunnel experiments on non- 
lifting cone-cylinder bodies in combination with wings of simple 
shapes were performed. Pressure distributions were measured at 
M = 3 and M = 4, both around the bodies and on the wings 
separately, as well as in combination, and it was found that the 
second-order interference was predicted reasonably well by the 
simplified theory. 


Symbols 
Cp = pressure coefficient, (p — po)/(1/2)pU? 
ACp = pressure coefficient in empty test section 
h = function defining the location of the (upper) 
wing surface 
M = free-stream Mach number 
p = local static pressure 
Po = free-stream static pressure 
R = hyperbolic radius 
U = free-stream velocity 
w = local (linearized ) upwash on wing 
i ae = Cartesian coordinates, the positive x axis in 


the free-stream direction 


¥ = ratio of specific heats 

65 = wedge semiangle 

A = leading-edge sweep 

p = free-stream density 

’= ¢+¢ = perturbation velocity potential 
¢ = first-order potential 
¢=wv+x = second-order potential 
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for) 


y = particular solution of nonhomogeneous wave 
equation 

x = complementary homogeneous solution 

Subscripts 

B = body 

z = interference 

m = measured values 

W = wing 


Introduction 


— SEARCH for vehicle designs which are aerody- 
namically efficient at supersonic speeds has aroused 
considerable interest in interference problems in recent 
years. In a large majority of papers on this subject 
linearized theory has been used. According to this, 
perturbation velocities and pressures due to interfering 
wings and bodies can be added linearly. The inter- 
ference thus affects the flow only through a change in 
tangency conditions. Linearized theory gives reason- 
ably accurate results for surface pressures for low super- 
sonic Mach numbers not too close to unity. At higher 
values of M, however, nonlinear effects become im- 
portant and it is no longer possible to add velocities 
and pressures linearly. Therefore, the wing-body flow 
field will be different from the sum of the separate 
wing and body fields, even in cases where tangency con- 
ditions are not directly affected. 

The linearized solution may be considered the first- 
order term in a series expansion, in terms of wing or body 
thickness ratio, of the full solution. The second-order 
term has been calculated in a number of cases (Buse- 
mann,' van Dyke?) and is known to improve the result 
considerably up to fairly large Mach numbers. Thus, 
for example, the second-order solution gives the pressure 
on a wedge of 5° semiangle correct within 5 per cent for 
M between 1.3 and 7, whereas the error of linearized 
theory is always larger than about 10 per cent, even 
for low M. 

In the present paper the problem of second-order 
wing-body interference at comparatively high Mach 
numbers (MM > 3) is considered. Since no practical 
second-order solution for the general three-dimen- 
sional case has yet been found, the problem is attacked 
by considering two simpler extreme models which are 
thought to bracket the true case. To check the theory 
a series of simple wing-body combinations have been 
tested at 7 = 3and 4. 
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Theory 


The problem considered is that of finding the pres- 
sure distribution on a wing situated in a disturbance 
field such as that created by a body of revolution or by 
another wing. To second order in perturbation ve- 
a perturbation velocity exists and is governed 
by the well-known second-order differential equation 
for the velocity potential @ due to small disturbances 
in a steady supersonic flow: 


locities, 


—(M* — 1)%,, + %, + $, = 
M%y — 1)%,(,, + &,, + &,,) + 
M2(0/dx) (6,2 + &,2 + 2) (1) 
(The undisturbed stream flows parallel to the positive 
x axis and is assumed to have a speed of unity.) From 
physical considerations, the boundary conditions of 
the problem require that the flow be tangent to all wing 
and body surfaces and that the flow perturbations 
vanish far upstream. Also, on shock surfaces the shock 

relations must be fulfilled. 

A familiar approach to the solution of this nonlinear 
equation for ® is an iteration procedure where ® is 
successively approximated by its first- and second- 


order parts. Setting 


P=o+¢ (2) 


where ¢ is the first-order, or linearized, approximation 
to @ and ¢ is the second-order term of the solution, the 
resulting differential equations are” 


—(M? — 1)grz + Gy + G22 = 0 (3) 
—(M on 1)@,. + f,, ig ®., _ 
M4(4 —_= 1) 9, ¢,77 + M?(0/O0x)(¢,? + $y? + $2") (4) 


In the following it is assumed that the linearized part of 
the solution is already known. For simplicity, we will 
only consider the case of a nonlifting body and wing so 
that the linearized solution may be written 
¢ = ¢p + ow (5) 

where ¢z is the solution for the body alone and gy that 
for the wing alone. For the second-order term we set 
¢ = 2 + bw t+ $i (6) 

where, similarly, ¢g and ¢w are the second-order terms 
for body alone and wing alone, respectively, and where 
¢; is the second-order term due to interference. Hence 
gp and ¢y are individually solutions of Eq. (4) with the 
appropriate subscripts for each right-hand term. 
Second-order solutions for plane and axisymmetric 
flow were considered by van Dyke.? Three-dimensional 
effects for thin wings at high Mach numbers were cal- 
culated approximately by Landahl.* Thus, further 
consideration of ¢g and @¢w is not required. After 
substituting Eqs. (5) and (6) in Eq. (4) and cancelling 
identities, the differential equation for ¢; is found to be 
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— (i? = 1) Pizz + Piyy = Piz: 


,of eh 
2M? Ox (1 + 9 ur) Sarews T 


To simplify the formulas we will from now on consider 
only such high Mach numbers that JJ? > 1, but at 
the same time let the body and wing thicknesses be- 
come small so that the hypersonic similarity parameter 
is low enough to justify the use of second-order theory. 
(For a simple wedge this means that the similarity 
parameter should not be higher than about 0.7, accord- 
ing to the Introduction.) The term ¢gw,¢g, will then 
be of order / tan A and the term ¢y-¢,- of order J/?. 
Restricting the analysis to wings of small sweep in com- 
parison with that of the Mach cone we may therefore 
neglect the former term in comparison with the latter, 
so that Eq. (7) simplifies for large 1/7 to 


YwySBy T HW wen: | 


(7) 


— M°$,,, T Piyy si Pizz = 


M?(0/0x)[M?(y — 1)¢ar¢wr + 2¢w-¢e:] (8) 


Eq. (8) is actually a useful approximation down to 
much lower values of 1/ than might appear from com- 
parison with Eq. (7) [note in particular that 1 has been 
neglected in comparison with (y — 1)47?/2.] It was 
shown in reference 3 that Eq. (8) leads to the piston- 
theory approximation, expanded to second order in 
perturbation velocities, and this is known to hold well 
for M > 3. 
Eq. (8) is the basis of the analyses that follow. 

high M/ the tangency condition for ¢, is simply that 


For 


G2 = —he,z-. forz = 0 (9) 


This arises from the tangency condition being applied 
on z = O instead of on z = h. All other second-order 
contributions to the tangency condition vanish for high 
M. 

The problem of finding practically useful solutions 
of Eq. (8) hinges on the possibility of discovering simple 
particular solutions in terms of gg and gy. For the 
general three-dimensional case such solutions have not 
yet been found. Instead, we will consider two simpli- 
fied models which are believed to approximate the true 
wing-body interference in two extreme cases. In 
both cases we will assume that the wing sweep is small 
in comparison with that of the Mach lines. Then, 
for high Mach numbers, ¢w can be approximated by 
the solution of 


— Mow + owe: = 0 (10) 
which is 
l diteng l 
~w=-r- at. : w(t, y)dé = — Vv h(x — Mz, y) 
(11) 


Here w(x, y) is the upwash required to fulfill the linear- 
ized tangency condition, gw. = w = h,, on the wing, 
and xz x. is the location of the leading edge at the span- 
wise station y. Eq. (11) is easily seen to be the linear- 
ized piston-theory approximation® of gy. 
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The simplifications introduced by Eqs. (10) and (11) 
are still not sufficient to make a simple solution of Eq. 
(8) possible. In the two flow models that are con- 
sidered we therefore make additional simplifying as- 
sumptions for ¢z. 

In the first model we assume that the curvature of the 
body field in the spanwise direction (y direction) is 
small so that gg is approximately a solution of 


— Mop, + ¢a:2 = 0 (12) 


This corresponds to a body field generated by two 
symmetrically placed wings as shown in Fig. 1. Two 
wings are required in order to correspond to the zero 
angle of attack case. Note that both the generating 
wings and the wing itself may be swept if the sweep is 
small in comparison with that of the Mach lines. 

Eq. (8) then simplifies to 


— Miz =] Pizz = 
M2%0/dx) [My — 


lL) ear¢wr = 2¢n:¢w:] (13) 
A particular solution, ¢; = y;, of this is easily found by 


use of results given by van Dyke? to be 


vi = M*[(3 — y)/4](0/0x) (gwen) + 
[M?(y + 1)/4]2(¢ar¢w: + ¢a:-¢wr) (14) 


To this must be added a homogeneous solution, x;, to 
satisfy Eq. (9). Using Eq. (11) we find that 


l x— Mz | 
xi = Vl f [Wi2(é, y, O) + h(é, V) Gp22(é, y, 0) \dé 
7 *L.E 
(15) 
The interference pressure may now be calculated 


from the second-order expression for the pressure 
coefficient which is” 

C, = —2¢6, + B°@,” _ ®,” it $,” (16) 
Introducing Eqs. (5) and (6) and considering only the 


second-order interference pressure on the wing at high 
Mach number gives 


2¢:r): 0 (17) 


ae = (2M? gp ews = 
Introducing ¢; = ¥; + x; from Eqs. (14) and (15) and 
¢w from Eq. (11) gives, after some simple calculations, 
Cot = —M(y 5 1) gp,W = [(y = = 1)/2|Mlgarz (18) 


For regions where gg, does not vary rapidly with x, 
or for regions near the wing leading edge, the second 
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Second approximate flow model. 


term in Eq. (18) may be neglected, leaving the very 
simple approximate formula 


Ges we —M(y = ia 1) oprW (19) 


It is easily shown that Eq. (19) is the second-order re- 
sult, at high V/, for the interference pressure that is ob- 
tained by applying the tangent-wedge approximation 
for calculation of the wing pressure, but using local 
values of pressure and Mach number (in the body field) 
instead of the free-stream ones. Thus, the present 
“local-flow tangent-wedge’’ method differs from the 
generalized shock-expansion method® which uses the 
values of pressure and Mach number at the leading 
edge. 

The approximations leading to Eqs. (18) and (19) 
may seem to make the results of little practical value 
since the body field in reality certainly seldom can be 
approximated as two-dimensional in the x, z plane. 
We therefore aJso will consider the other extreme 
namely, when the body field is approximately two- 
dimensional in the x, y plane, as illustrated in Fig. 2, 
that ‘s, when the curvature of the body field normal to 
the wing plane may be neglected. For this second 
approximate flow model we should thus solve Eq. (8) 
with vn, = 0; i.e., 


— Moir, + Piy, + Pizz at 
M*(y — 1)(0/0x)(¢er¢w:) (20) 


It can easily be verified that a particular solution of this 


equation is 
Vi = —(] '2)M?(¥ - 1)(O Ox) (¢Gwe¢p) (21) 


The homogeneous solution x; should in this case satisfy 
the three-dimensional wave equation 


— MX zz = Xin + Xizz = 0 (22) 


and should cancel the normal derivative of y; on z = 0; 
1:€;; 


Xiz2 = (1/2)M?(y — 1)(0/0x)(¢sew:) onz =0 (23) 


It is well known that a solution with these properties 


is given by 
xi = — (1/m) SS (1/R)xilé, ndé dn (24) 
where 


R= V(x — 8? — M[(y — 9)? +27] (25) 
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and where the surface integral is to be extended over 
the region of the é, 7 plane intercepted by the forward 
Mach cone from the field point (x, y, 2). 

However, 


$=Vvit xi (26) 


thus determined may not necessarily make the shock 
relations on the oblique shocks from the body and wing 
fulfilled to second order. It is possible to simplify these 
relations by applying van Dyke’s smoothing proce- 
In this, shocks are considered as the limit of 
This can be 


dure.” 
infinitely rapid isentropic compressions. 
achieved by assuming the wing and the two-dimen- 
After the 
analysis is carried out these cusps are shrunk to zero. 


sional body to have cusped leading edges. 


Let us assume that the body field, after the smooth- 
ing procedure has been applied, is such that ¢, and its 
first-order derivatives vanish for x < My. Eqs. (21), 


(24), and (26) will then give 
dic = (1/2)M?(v — lewearr(0, 0) onx = My (27) 


i.e., the interference pressure will not vanish for a 
vanishingly small body-pressure perturbation. How- 
ever, the revised characteristic in the wing field is found 


to be given by 


x= My + [(y = 1) 2|\Meow (28) 





and 


iG: 


(3/2 —f)r 


F(t) = 
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in the present high-/ approximation. On the surface 


defined by Eq. (28) 


YBr = Pir = (x — My) ¢erz(0) = 
(1/2)M?(y — l)gwearr(0) (29) 


is found to vanish to second order (this may require 
analytic continuation of y,, for x < My) so that on the 
true characteristics from the “‘body’’ leading edge the 
infinitesimal pressure perturbation from the body will 
give a zero total pressure perturbation as required. 
Similarly, one can show that on the true characteristics 
from the cusped wing leading edge the total pressure 
perturbation will vanish. Hence, the solution given 
above is complete. 

We consider first the following case, where the body 
field accelerates linearly from x = My with ¢z,, 
Cp, and the wing is a simple wedge—namely, 


) 


gp = (1/2)Ca(x — My)? (30) 
gw = (6/M)(x — Ms) (31) 


where 6 is the wedge semiangle. For the smoothing, an 
infinitesimal cusp must be added to the leading edge 
of the wedge, but this will only have influence in Eq. 
(24). After some calculations, the equations listed 
above give for the interference pressure on the wing 

Cyi = —2M Cpx-6 — [M(y — 1)/r]Cpx-6- F(t) (32) 
where t = My/x 


fort > — 1 
fort < — 1 


Since the interference pressure is linear in ¢, this result can be used to build up the solution for an arbitrary body 


field by superimposing the interference pressure due to infinitesimal body fields defined by 


d 


_ f(1/2)(x — & — My)? denzr(E, y) 


forx —&> My 


_ L0 forx —§< My 
After an integration by parts the result may then be expressed as follows: 
Cy = —M(y + 1)¢86 + M(y — ve f gar(&, y)G(x, &)dé (34) 
where G(x, ) = — (1/m)[e — cos—! (é/x) — (2 — t/x)V (x + 8/(x — 8] (35) 


The limits of integration follow from the fact that the interference pressure is influenced only by those body perturba- 
tions over the wing which are inside the upstream Mach cone from the field point, and from the fact that the body 


field is a function of x — My only. 
Stiltjes integral. 


The smoothing means that at jumps in ¢z,, Eq. (34) should be interpreted as a 


Finally, by superimposing results for infinitesimal wedges whose first-order flow fields are defined by 


dow = 10 


the interference pressure on an arbitrary wing may be 
obtained as follows: 


Coro = —M(y + 1)¢n.~ + M(y — 1) , Wz,(x)dx, X 


if on pee (E)G(x —%x%,& — wa | (37) 


which also should be interpreted as a Stiltjes integral 
at jumps in w and gz,. 


If gg, varies very little with 


— J-CU/M)(x — x, — Mz)dw(x) 


forx — x, > Mz 


36 
forx —x,< Mz aie 





x the simplified result Eq. (19) is again retrieved. 
Comparisons between the different formulas are 
given in Fig. 10, together with experimental results. 


Second-Order Wind-Tunnel Correction 


Assume that calibration of a supersonic wind tunnel 
has given as a result a pressure coefficient distribution 
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in the empty test section, AC,, referred to some average 
or nominal Mach number Jf. A complete calibration 
would also require measurements of local flow inclina- 
tions, but for the present purpose we shall assume that 
these are negligible. The linearized, or first-order, 
correction to the measured pressure distribution is 
simply 

C, = Com — AC, (38) 


where C, stands for the corrected value and C,,, for 
the uncorrected one. Using Eq. (19) we may derive a 
simple second-order correction formula by assuming 
that the flow nonuniformities in the wind tunnel com- 
prise the “body” field. To first order we may then 
replace gg, in the formula by —(1/2)AC, and, for a 
wing-like member, we may set w = (1/2) MC,. Hence, 
applying Eq. (19) gives 


Com = C, + AC, + (1/4)(v + 1I)M?AC,C, (39) 
or, solving for C,, 


. Cun — AC, ; 
C, = ) SIRES (40) 
1 + (1/4)(y + 1I)MPAC, 

Alternatively, one can use local values of static and 

dynamic pressures in the empty test section when com- 

puting the pressure coefficient, instead of average values 

as above. By some simple calculations one can show 

that for small velocity errors at high Mach numbers 

this procedure corresponds to replacing Eq. (40) by the 
formula 

¢, « = Se. (41) 

1 + (1/2)M?AC, 

The difference between Eqs. (+1) and (40) is generally 

small. Eq. (40) was used for correcting all measured 

results presented in the present paper. At 1J = 4 
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the second-order correction amounted in some cases 
to 15 per cent of the measured pressure coefficient. 


The Experimental Investigation 


The experimental investigation was performed in 
two small wind tunnels of the vacuum type, with 
atmospheric stagnation conditions. The pump capac- 
ity for the vacuum chamber—which is mainly used 
for much larger tunnels—is sufficient to run the small 
wind tunnels continuously. The Mach 3 wind tunnel 
with a test section of 24 X 30 cm.” has a nozzle giving 
a mean value of Mach 3.02 with a rather good velocity 
distribution. The Mach 4 wind tunnel has a test sec- 
tion of 16 XK 20 cm.*. The mean Mach number ob- 
tained in the test section was 4.01, but the velocity dis- 
tribution cannot be said to be quite satisfactory, since 
the local Mach number varied between 3.95 and 4.15. 
However, the formula worked out above [Eq. (40)] to 
correct for nonuniformities in the test section seems 
quite capable of dealing with the distribution obtained. 
The influence of the nonuniform wind-tunnel flow is 
demonstrated in Fig. 16. 

The wing models used consisted of three delta wings 
and one straight wing, which all could be attached to a 
cone-cylinder body at different positions aft of the 
shoulder, thus being located in a nonuniform flow field, 
including both compression and rapid-expansion re- 
gions. The strength of this field could be varied by 
using different cone-angles for the conical noses. Three 
different cones, with cone semiangles of 15°, 17.5°, 
and 20°, were available. The straight wing had a 
wedge profile, the half-wedge angle being 4.1°. Its 
chord length and span are given in Fig. 3, where the 
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plan-form parameters for the delta wings also can be 
found. The leading-edge sweep angles were chosen 
to be 30°, 40°, and 50°. The span was the same for 
all wings and equal to three times the radius of the 
cylinder. The profile for the rear edge section of the 
delta wings consisted of two parabolic arcs, the vertex 
located at the root chord. The maximum thickness for 
all wings was located at the trailing edge and had for 
the symmetry sections the same absolute value for the 
four wings. This value was taken to be 10 per cent 
of the root chord for the 50° swept wing. The surfaces 
of the delta wings were generated by straight generators 
from the apex. 

The four wing models, the length of the cylindrical 
part, and the conical noses were easy to change accord- 
ing to the design demonstrated in Fig. 4. 

The wings had pressure holes located on one surface 
along four sections parallel to the axis of symmetry. 
For two representatively located pressure holes there 
were also holes on the opposite surface in order to 
check the zero-angle setting. 

To get the basic pressure distributions for the four 
wing models in uniform flow a reflection plate was used. 
This plate shielded the flow field around the wings 
from the nonuniform flow created by the conical noses. 
The dimensions of this plate and its position are shown 
in Fig. 3. 

For convenience, the straight wing is hereafter de- 
noted by W, and the delta wings corresponding to 
30°, 40°, and 50° sweep angles by We, W3, and W,, 
respectively. 

The static pressures were recorded on a mercury 
multitube manometer for the M = 3 investigation 
and on a silicone-oil manometer for M = 4. Pitot 
pressures were always recorded on the mercury- 
manometer. 
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The experimental investigation included, in addi- 
tion to static pressure distribution measurements for 
the different cone-cylinder-wing configurations, cali- 
bration of the sections by means of pitot tube measure- 
ments; this also was the method of determining the 
flow field above the reflection plate and around the 
cone-cylinder. In the latter case the results were 
evaluated under the assumption that the pitot meas- 
urements corresponded to zero flow inclination between 
the streamlines and the pitot tube, the latter always 
being aligned along the axis of the symmetry. It had 
been confirmed that the pitot tube was insensitive to 
variations in the alignment of the flow within limits 
not exceeded in the present investigation. 


Results and Comparisons With Theory 


For brevity, only some results for the 15° cone- 
cylinder and the wings W;, W2, and W; are presented. 
Results from the body-field pressure measurements 
for the 15° cone-cylinder are shown in Figs. 6 and 7. 
At M = 3, theoretical results were obtained by the 
method of characteristics. As seen, the agreement be- 
tween theory and experiments is excellent, except in 
Sections 2 and 3 far behind the conical shock. This 
small discrepancy is probably due to poor numerical 
accuracy in the theoretical results for regions influenced 
by many Mach line reflections on the cylindrical part 
of the body. The solid curves shown in Fig. 7 for M = 
4 were obtained by applying the unified hypersonic- 
supersonic similarity law‘ to the mean experimental 
results for the 20° cone-cylinder at M = 3. The 
favorable comparison in the whole expansion region 
provides an interesting check of the similarity law. 

Wing-alone results for the wings We and Ws; are 
shown in Figs. 8 and 9. The theoretical curves were 
obtained from Lighthill’s third-order piston theory° 
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with higher-order terms in 1/7? (including three- 
dimensional effects) calculated from Landahl’s theory.* 
It was found that the boundary-layer effect on the pres- 
sure distribution was rather large, so that when cal- 
culating the wing pressures (as well as when calculating 
the second-order interference) the wing thickness was 
increased by the displacement thickness of a flat-plate 
compressible laminar boundary layer. (The Reynolds 
number, based on local wing chord, ranged from 0.5 X 
10° to 5 X& 10°.) As seen, the agreement between 
theory and experiments is quite good when the bound- 
ary-layer thickness effect is included, except for Section 
1, which is partly located inside the region bounded by 
the Mach line from the wing apex where the theory 
used does not apply. The measured pressures at the 
holes closest to the trailing edge are consistently some- 
what lower than the theoretical values. The reason 
for this is not quite clear, but these holes are all located 
so close to the trailing edge (0.6 to 3 mm.) that the 
pressures might be influenced upstream through the 
boundary layer by the low pressure at the blunt 
trailing edge. 

For the pressure measurements on the wings in the 
body field, results are presented only for those cases 
for which the interference was found to be large. 
Theoretical results at J = 3, obtained using Eqs. (18), 
(19), and (34), are shown in Fig. 10 together with the 
measured values, for Section 3 of wing W, in its forward 
position. In evaluating the theoretical results the 
linearized approximation gz, = — (1/2)C,g was used 
throughout. As may be seen, the different approxima- 
tive results differ very little from each other in view of 
the large body-pressure gradients present. Also, the 
simplified formula Eq. (19) gives results which lie in be- 
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tween those of the two extreme flow models. As per- 
haps would be expected, the second flow model gives 
the best agreement with the measured pressures. The 
reason for the large discrepancy at the point closest to 
the leading edge is not yet known, but it is suspected 
that shock-wave boundary-layer interaction at the in- 
tersection of the body shock with the leading edge may 
be involved. In any case, Fig. 10 suggests that the 
effect of body-field velocity gradients in general is small 
and that the proposed local-flow tangent-wedge approxi- 
mation, given to second order by Eq. (19), should give 
small errors at high Mach numbers, even for regions 
with large body-field gradients. 

Results for Section 4 of the wings We and W; at 
M = 3, in positions as indicated by the sketches, are 
shown in Figs. 11 and 12. Apparently the agreement 
between theory and experiment is good near the leading 
edge but the experimental values consistently decrease 
more rapidly towards the trailing edge than the theory 
predicts. 

Results for WM = 4 are given in Figs. 13-16. The 
agreement between theory and experiment for this 
Mach number is everywhere excellent, except on the 
rear portion of Section 4 where the theory overestimates 
the interference. A satisfactory explanation for this 
discrepancy has not yet been found. It follows from 
Figs. 13 and 14 that the interference pressures are of 
the same order of magnitude as the body pressures. 
This shows that the linearized theory. can seriously 
underestimate the wing-body interference at high Mach 
numbers. 

In Fig. 15 a comparison is made between results 
obtained from the generalized shock-expansion method,® 
the present local-flow tangent-wedge method, and ex- 


periments. Apparently the Mach number is not 
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sufficiently high (or the body field varies too rapidly) 
for the shock-expansion method to be applicable in this 
case. 

Fig. 16 is prepared to illustrate the influence of the 
second-order wind-tunnel correction. In order to 
demonstrate the importance of the second-order effect, 
that section was chosen (Section 3) which had the larg- 
est wind-tunnel correction. As seen, the second-order 
correction is everywhere larger than the first-order 
remarkable improvement 


correction and causes a 


between experiments and theory. 


Conclusions 


By neglecting the curvature of the body flow-field 
in the spanwise direction, or in the direction normal to 
the wing surface, two approximations were obtained 
for the second-order wing-body interference at high 7 
on wings with small sweep in comparison to that of 
the Mach line. Despite the entirely different under- 
lying assumptions, the two approximations were found 
to give results that differed very little, even in regions 
of rapidly varying body field, and for most practical 
cases they can probably be replaced by the very 
simple approximate formula, Eq. (19). This corre- 
sponds to using a method for calculating the pressures 
on the wing that might be termed “‘local-flow tangent- 
wedge” approximation and which is identical to the 
ordinary tangent-wedge method, except that local 
values of pressure and Mach number in the nonuniform 
body field are used instead of free-stream values. The 
experimental investigation confirmed that this simple 
method predicted the second-order interference well 
in most cases. At \/ = 3 the agreement between theory 
and experiments was reasonable and at MJ = 4 gener- 
ally excellent. The less good agreement at MW = 38 is 
probably due to the combined effect of low Mach 
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number, wing sweep, and a more rapidly varying body 
field than at Jf = 4. A comparison with the general- 
ized shock-expansion method* at .1J = 4 shows that 
the latter theory gives results that are in poor agree- 
ment with experiment and with the present theory, 
except very near the leading edge. Apparently, the 
Mach number is too low for this method to be appli- 
cable. 

To second order, the tangent-wedge method is, of 
course, identical to the shock-expansion method or to 
Busemann’s formula. Hence, the present results indi- 
cate that the pressure at a particular point on the wing 
in the nonuniform body flow is very nearly the same 
as if the wing were immersed in a uniform stream with 
flow properties equal to those of the body flow at the 
point in question. It would be worth while to investi- 
gate whether this simple result also holds for wings 
with blunt leading edges and for blunt-nosed bodies, 
for which there appear large entropy gradients in the 


flow field. 
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Multicondition Terminal Control Svstems 
for Aircraft’ 


RICHARD V. MORRIS* anp CARL W. STEEG, JR.** 
Massachusetts Institute of Technology 


Summary 


A method is presented for synthesizing a type of control sys- 
tem applicable to automatic landing of aircraft and designated 


multicondition terminal control system. This closed-loop 


is 
control system is designed to cause the response of a controlled 
member to meet several independently specified conditions at a 


particular instant of time. The desired values of altitude and 


altitude rate at touchdown typify such conditions. Because 


the system operates to shape the response at one time only, it is 


1€ 


suitable for use with controlled members that have oscillatory 


or divergent responses to forcing inputs 


Symbols 


matrix of coefficients a;;(t) relating response vari 
ables and their derivatives 
coefficient relating 7; and x; 


matrix of coefficients 5;(t) relating forcing to 


B = 
derivatives of the response variables 
t = coefficient relating 7; and y 
C1, ¢ = parameters of linear feedback element in follow 
up servo 
C3, ¢ = parameters of linear controller element in follow- 
up servo 
T2T = weighting function relating x;( 72) to x;(7)) 
H(t = weighting function of superposition integral used 


in determination of x;(72) 
H = jnitial aircraft altitude, ft. 


hj(72, 7 = weighting function of superposition integral used 
in determination of x;(72) 

I T,t = integral of y weighted by /; between ¢ and 7 

1,J = dummy variables of summation 

Ky = gain of computer-study aircraft model 

Ky = gain of linear controller element in follow-up 
servo 

m = one less than number of conditions on response 
components 

n = order of set of differential equations 

t = time, sec 

¢" = time at which xg, becomes zero and assumed 
forcing program is applied 

7 = touchdown time, end point of the control interval 
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Tc(A = transfer function of linear controller element im 
follow-up servo 

Tr(A = transfer function of linear feedback element in 
follow-up servo 

l = aircraft forward speed, ft./sec 

>. 4 = response vector with x; as components 

W: = actuating error 

Cs, %y = response components, components of X 

xia(T,t = estimate of value of x;(7) made at ¢ 

v;D = desired value of x,;(T 

v;DI = predicted value of x;p 

\ = value of x; att = 0 

X26 = error in x2(7) caused by input limiting 

b 4 = forcing vector, with y as nonzero component 

= forcing input 
‘ = assumed forcing program 

VI = limiting amplitude of y 

r = aircraft flight-path angle, rad 

c = damping ratio of computer-study aircraft model 

r = complex frequency variable or Laplace transform 
variable 

o = time constant of computer-study aircraft model, 
seconds 

T = dummy variable of integration 

Ti, 72 = particular times 

> = fundamental solution matrix to general vector 
differential equations 

? = limited assumed forcing program 

a = computed value of assumed forcing program 

oc; = adjustable constants in assumed forcing program 

PCN = numerator in fraction determining ¢¢ 

y¥ = known function of time in assumed program 

vy; = known functions of time in assumed program 

w = natural frequency of computer-study aircraft 
model, rad. /sec 

sgn = algebraic sign 


Dot denotes differentiation with respect to ¢ 


Introduction 


F  ieninonse TOUCHDOWN is an example of a motion- 
control problem that frequently confronts aero- 
nautical engineers: A body is to be moved in such a 
way that both its position and velocity coincide with 
desired values at a known instant of time. As its 
forward velocity carries the aircraft across the touch- 
down zone on the runway, both lateral and vertical 


NPuUT ine AC ACTUATING FORCING vstew 
ata COMMANE NPUT ERROR NeuT NTROLLE RESPON 
a | ° a i nal NTROLLER ~ * . 
ExTRA s + MEMBER 
. 
MMANDED 
RESPONSE 
L FEEDBACK le. 
MEMBER 
Fic. 1. Block diagram of a basic automatic control system. 
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position and velocity are to assume the desired values. 
The lateral position is to coincide with the center of 
the runway at a zero rate of change. The altitudes of 
the aircraft and pavement are to coincide at a small 
sinking speed for the aircraft. 

In the past, most touchdown control systems have 
been designed to cause the aircraft to meet these condi- 
tions by attempting to follow a specified final-approach 
path. To ensure that the touchdown conditions are 
met, these systems overspecify the motion during the 
final approach. As a result of this overspecification, 
attempts to minimize path errors that occur before 
touchdown may prevent attainment of touchdown con- 
ditions. The subject of this paper is a method of 
automatic motion control that stresses the importance 
of position and velocity at touchdown by eliminating 
the requirement of continuously minimizing a path 
error. The method is designated multicondition 
terminal control. The principles of operation of this 
method of control and the use of a mathematical state- 
ment of the principles to aid in the synthesis of con- 
trol equipment are discussed in following sections of the 
paper. Analog computer results are used to demon- 
strate the performance of a physical control system 
synthesized in the manner discussed. Random-signal 
disturbances and other variations from design condi- 
tions are introduced in the computer tests 


General Description of Terminal Control 


Terminal control systems and conventional servo- 
mechanisms comprise two distinct classes of control 
systems. The present section relates terminal control 
to more familiar methods in terms of some basic con- 
cepts of automatic control. A general block diagram 
for an automatic control system is shown in Fig. 1. 


The distinguishing characteristic of 
automatic control systems is the presence of an actu- 
ating error that controls the system response. A 
forcing input to the controlled member is generated 
from the actuating error by the controller. To ac- 
complish its purpose, an automatic control system must 
be supplied with data containing information pertinent 
to the response desired by the user. The operation of 
command extraction provides this information in an 
appropriate form as the command input. The actu- 
ating error results from a comparison of this input with 
the command response, a function of the system re- 
sponse generated by the feedback member. 


closed-loop 


Performance specifications are a necessary prerequi- 
site to the determination of control system configura- 
tions that will cause the motion of the controlled body 
to meet the user’s requirements. The specifications 


are based on a detailed description of the motion of 
the body, but often more than one such description 
defines motion that will meet the user’s requirements. 
In such instances, the first task of the system designer 
may be the study of the various descriptions to deter- 
mine those useful as a basis for realistic performance 
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specifications. A quantity that aids in the formula- 
tion of the specifications is the response error, defined 
as the difference between the desired and actual re- 
sponses at the same instant of time. 

For example, the primary task of a touchdown con- 
trol system is to cause the aircraft to meet the touch- 
down conditions, which can be assumed to occur at a 
known time after the control system is actuated. Em- 
phasis is placed on minimizing the response error at 
the touchdown time 7. The control system can be de- 
signed to weight the value of the response error at 7 
more heavily than previously occurring values. A 
performance specification that might be used is that 
the response error at 7 be the only value considered, 
and all efforts are bent on minimizing this value. Such 
is the action of a terminal control system. 

Time-varying operations are necessary in the com- 
mand extraction and feedback members because the 
required estimations are made over a continuously de- 
creasing interval. Techniques have been developed 
for deriving configuration of these members so that 
only linear operations are performed on the input sig- 
nals. The techniques are reviewed in Appendix I. 
Although the fact may not be obvious from the intro- 
ductory discussion, simultaneous control of the posi- 
tion and velocity at touchdown is made possible by 
proper choices of the system configuration and the as- 
sumption relating to the future forcing. Because the 
control system operates to shape the system response 
at only a particular instant of time, it is suitable for 
use with controlled elements that have oscillatory or 
divergent responses to a forcing input. 


Type of Controlled Member Necessary 


The method of automatic control discussed here is 
applicable to controlled elements that are linear and 
have a single control input, which may be limited in 
Several formulations of the mathematical 
The one chosen 


amplitude. 
model of such a member are possible. 
for this paper is a set of m simultaneous linear first-order 
differential equations. Such a model represents any 
constant-coefficient linear system describable by a 
transfer function that is a rational fraction in \, the 
complex frequency variable. In addition, the model 
represents time-varying linear elements describable 
by a single differential operation of order m on input 
and response, provided that the order of the part of the 


operation on the input is less than 7. 
n 
= > ay(t)x; + b(t)y, «= 1,2,...,” (1) 
j=) 

Furthermore, the variables x; usually can be related 
more closely to measurable variables in the problem 
than can the high-order derivatives that are involved 
in a single operation of order m derived from the set of 
Eq. (1). The set of equations is adaptable to matrix 
manipulations that can be used to derive several useful 
results. 
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MULTICONDITION 


Principles of Operation 


The principles of operation of a terminal control 
system for simultaneous control of a response and its 
derivative are stated in nonmathematical terms in the 
present section. The following section uses a mathe- 
matical statement of the five principles of operation as 
a basis for the synthesis of control equipment. 

A fundamental difference between the terminal con- 
trol and follow-up servo systems is the nature of the 
actuating error. In a follow-up system this error is 
approximately the present value of the response error. 
In a terminal control system the actuating error is the 
difference between the estimated values of the actual 
and desired system responses at 7, an estimate of the 
final response error. The determination of the esti- 
mates is one phase of the problem of linear prediction.': * 
For the estimation of a quantity, the input information 
to the estimation device may be regarded as the result 
of applying some member of a class of inputs to a shap- 
ing element with known characteristics. To construct 
an effective estimation device, knowledge of the char- 
acteristics of the shaping element and an assumption 
relating to the class of inputs to this element are neces- 
sary. An example is the assumption in the applica- 
tion® of Wiener’s work that the signals used by the 
estimation device are members of a particular statistical 


ensemble. 


Synthesis of Components 


The five principles of operation are examined here 
individually to determine a synthesis procedure. Cer- 
tain mathematical results, which are given without 
Details of the derivations are 


proof, are interpreted. 
Because certain components 


given in the appendixes. 
of the control system depend on the form of the con- 
trolled member, the following simple mathematical 
model, a special case of Eq. (1), is used to illustrate the 
discussion : 


x1 = Xo \ ( 


Xx. = 


2) 


—x2 + yf 


The general solution to a set of linear differential 
equations provides an adequate mathematical state- 
ment for the synthesis of the estimation devices. This 
solution is derived in Appendix I; the result is Eq. (3). 
The details of the derivation® are not necessary for the 


present discussion. 


n 72 
Xi(7T2) } > fij(t2.71 )xj(71) + f hj( tT, T) V(r) dr 
j=l Tl 


+= 1,2,...,2 (8) 


The equation states that x;(72), the value of a particular 
response at time 72, is determined completely by the 
values of all the system responses at another time 7 
and by the forcing applied to the system between 7; 
and 7. The relationship between x;(72) and the sys- 
tem responses at 7; is determined by the homogeneous 
solutions to the differential equations. The relation- 
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ship between x,(72) and the forcing is expressed as a 
weighted integral, or convolution, of the forcing with a 
weighting function related to the homogeneous solu- 
tions. The functions /;;(72,7;) are related to the 
solutions to the homogeneous portion of the simul- 
taneous set of Eq. (1). These functions have the prop 


erty 


I ¢3( T2,T1) |eome Jl, : “J (4) 

rit (0, @#j 
A more detailed interpretation of Eq. (3) is that the 
variables x; interact, as specified by the coefficients of 
Eq. (1). Because the system is linear, the behavior 
of each variable x 
several terms, each involving only a single response 
One term is the behavior that 


can be considered as the sum of 


variable or the forcing. 
x; would have if it were the only nonzero response 
variable at time 7; and the system were unforced. 
This contribution is expressed by x;(11)fii(72,71). 
Another set of terms, x;(71) fi;(72, 71) with 7 ¥ 1, ex- 
presses the behavior of x; that would occur if each of 
the other x; were the only nonzero response at time 7; 
and the system were unforced. The final term in 
Eq. (3) is the superposition integral that represents 
the change in x; due to the application of the forcing 
input y during the interval between 7; and 7». 
For the example of Eq. (2), 


X1(71) +(1- os os X2(71) + 


{ (1 — e” ™)y(r) dr i 
Tl ; (2) 
e'  x0(71) + { Zs 


We now discuss synthesis of the portion of the con- 
trol system determined by Principle 1, ‘““‘The control 
system continuously estimates the final values of the 
These estimates 


X1( T2) = 


" y(r) dr 


Xo( T2) 


system response and its derivative. 
are based on a knowledge of the dynamic characteris- 
tics of the controlled member and on an assumption 
concerning the program of forcing expected over the 
remainder of the control interval.’’ This principle 
requires an assumption relating to the forcing to be 
applied between the present and terminal times, ¢ and 
T, respectively. The assumption is that the forcing 
is described by 


valt) = bc¥(7) (6) 


The functions f;;(7,t) and h,(7,r) are sufficient to de- 
scribe the controlled element for the purposes of con- 
structing the estimation device. 

The mathematical expressions for the estimates valid 
at the present time can be obtained from Eq. (3) by 
assigning the values 

m1=t a 1 { (7) 


A = Ja xX ;(T2) -_ xi4(T,7)§ 


When these values are substituted into Eq. (3), the 
estimated responses are 
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Fic. 2. Generation of the estimation signal x. 


n 


T 
Xia(T,b) = >» fis(T,D) x(t) +f h(T,r)ya(r)dr (S) 
t 


j=l 


In particular, for the controlled element used as an 
example, 


x4(7,t) = (fH) — ee xe(t) + 
7 
{ (1 — e”~") y4(r) dr\ (9) 
. 
Xeg(I,t) = * ati xo(t) + y? yvalr) dr 
t 


It is convenient to assume that electrical signals 
representing the values of x;4 are desired and that elec- 
trical signals representing the x,(¢) are available. For 
the computation of an x;, signal, each of the x, signals 
is acted on by the appropriate time-varying gain, 
fij(7,t), and summed with another signal representing 
the superposition integral, as indicated in Fig. 2 for the 
estimate x;4 of the example. The method of gener- 
ating the signals representing superposition integrals 
is deferred until Principle 5 is discussed. 

With the generation of the estimation signals estab- 
lished, the synthesis procedure can progress to include 
Principle 2, ‘The control system compares the desired 
and estimated final values of the system response to 
obtain the actuating error.” The comparison is made 
with a summing circuit, as shown in Fig. 3. In some 
cases, such as carrier landings, exact determination of 
the desired terminal values of the response components 
is not possible before ¢ = 7. In such a situation, an 
estimate of the desired terminal values x;pp must be 
substituted for the actual desired value x;p. The esti- 
mate is obtained from a predictor of some nature 
utilizing as input data the information available con- 
cerning X;p. The design of this predictor is a problem 
separate from the control 
systems.! ? 

The difference between estimates of the actual and 
desired values of the response components at ¢ = TJ is 
obtained in the same manner for each component. 
Because the terminal control system under discussion 
is designed to control two response components, a con- 
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venient assumption is that these two components are 
designated x; and x2, respectively. On this basis the 


single actuating error is 


Xp(t) = Xipp(t) ated Xia (1 ,t) 10) 


) 


Principle 3 defines the operation of the controller for a 
nonzero actuating error: ‘“‘If the actuating error is not 
zero, the control system operates on this error to gener- 
ate a forcing input for the controlled member that will 
reduce the magnitude of this error.’’ The 


change of magnitude of the error must be negative. 


rate of 


Appendix II demonstrates that this condition is equiva- 
lent to the condition that the rate of change of the 
error be opposite in sign to the error itself. Appendix 
II also shows that the largest negative rate of change 
of the magnitude is obtained if the conditions of Eq. 
(11) are met. 


l-y(t)| = |ya(t) 
sgn y(t) = [sgn /,(7,t)| [sgn xz(t)] 


( | ] ) 
lv(t)| — |y,(t)| is maximized 


No controlled member is completely linear, although 
a linear description is valid over a portion of the oper- 
ating range. In the nonlinear range the controlled 
member frequently saturates whenever the magnitude 
of y exceeds a maximum value y,. Because saturation 
is one of the principal nonlinear effects limiting the per- 
formance of controlled members, this representation is 
useful for many actual systems, including the ones dis- 
cussed here. Amplitude limiting also may be imposed 
on the control input, in an attempt to prevent undesir- 
able response effects such as aerodynamic stall. For 
controlled members with an input limit, an assumed 
program with an amplitude larger than the limiting 
value y, is not permitted because this forcing invali- 
dates the linear description of the controlled member. 
The first and third conditions of Eq. (11) are met if 


y(t) = 
yValt) 


(12) 


The second condition of Eq. (11) states that for a non- 
zero error the sign of y must be positive if the signs of 
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Fic. 3. Determination of the actuating error. 
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the actuating error and the weighting function /,(77)f) 
agree, and negative if they differ. Hence, the forcing 
input is determined partially by the switching opera- 


tion specified by 


a fty.,, sgn y(7\t) = sgnxg, Xe ¥ 0 (13 
Vl 5 or ; ”) 
; \—y,, sgn hy(7\t) ¥ sgnxg, Xe # O 


The remaining characteristic of y is determined by 
Principle 4, ‘‘When the actuating error becomes zero, 
the control system generates the assumed forcing pro- 


gram and applies it to the controlled member.” This 
principle requires that 
y(t) = valb) (14) 
if 
vz = 0 (15) 


The switching element to determine the forcing in- 
put can be a three-position switch controlled by signals 
representing the actuating error x,(¢) and the weighting 
function /,(7\t), as shown in Fig. 4 for the example 
system. The mathematical specification of the switch- 
ing operation is 


+ yp, sgn xg = sgn (1 — e*) 
vi Yall), Xe = O, (16) 
=¥ 
Fi — yz, sgnxg ¥ sgn (1 — e ) 


In practice, a small dead zone may be added to the 
switch to reduce chatter caused by noise in the control 
system. Then the forcing will be the assumed pro- 
gram if the magnitude of the actuating error is less 
than the dead-zone limit. 

The assumed forcing program is determined by Prin- 
ciple 5, ‘““The assumed forcing program is so chosen 
that the desired and actual values of the system response 
rate coincide.’’ The assumption made concerning the 
forcing program in Eq. (6) is helpful in synthesizing 
a device to meet the requirement of this principle. 


T 
f h<T,r) W(r) dr 
t 


is useful for the discussion. Substitution of the re- 
sults of Eqs. (6) and (17) into Eq. (8) yields 


The definition 


TAT SD) = (17) 
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xia(T\t) = D> fil Tit) x(t) + bel (Tit (18) 
j=l 
To satisfy the requirement, it is necessary to make 
Xe4(1 ,t) Vep (19 
, ») 


After substitution of Eq. (19) into Eq. (18) for 1 2, 
the resulting equation is solved to determine 


Xep — yo fo(7\t) x,(t) 


i=] 


I.(7,t 


(?0) 


Pc 


All of the terms on the right-hand side of Eq. (20) are 
known at the present time /, and ¢¢ is computed con- 
tinuously from presently available information. For 
the example of Eq. (2), 


T.(T,t) 


(21) 


A method of computing ¢¢ is illustrated in Fig. 5 for the 
example system. The difference between x2» and x2, 
is determined, and this difference is multiplied by the 
reciprocal of /.(7,t) to obtain ¢¢. Because /2(7,t) 
may become zero at certain instants during the control 
interval if the controlled member is oscillatory, an 
approximation to the reciprocal of /, that remains finite 
must be used. Such approximations usually have little 
effect on the system performance. A 
situation arises when Eq. (20) specifies a magnitude of 
oc larger than y,. Then, the estimate of x,,4(7,t) 
specified by Eq. (18) is not correct, although the actu- 


more serious 


ating error may become zero because the forcing applied 
to the controlled member is limited in magnitude and is 
different from the value of forcing assumed for the esti- 


mation. The difficulty can be lessened if 


V(r) = 1 (22) 


because the assumed forcing program is a constant of 
value ¢@c. If this signal is passed through a unity-gain 
limiter set to limit at y,, a new assumed forcing pro- 


gram results that is defined by 


(yz, Yt < bc 
74) = Qc, ak < de < VL (23) 
ly, bc0< Me 


If @ replaces ¢¢ in the calculation of the estimated re 


Fic. 5. Computation and use of the assumed forcing program 
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sponse x;4(7'¢), then the response condition will be met 
whenever possible, even if a value of |¢¢| > yy indi- 
cates that the rate condition cannot be met. In such 
cases, the rate error will be 


Xe = (Gc — O) IT t’) (24) 


where ¢’ is the time at which x, becomes zero and as- 
sumed forcing program is applied. 

If dc exceeds the limit level or if the actuating error 
does not become zero during the control interval, the 
terminal conditions are incompatible with the dynamic 
-apabilities of the controlled member when the type 
of forcing generated by the control system under dis- 
cussion is used. 


Experimental Results 


The controlled member used in the analog computer 
tests is a rough model of the vertical motion of a land- 
ing aircraft. If the aircraft is assumed to be trimmed 
to fly along a straight flight path inclined to the hori- 
zontal at an angle I’, the altitude response is approxi- 
mately 


25) 


H(t) = Ay + Uot tan [ 4 x1 
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and its derivative, the sinking speed, is approximately 





H(t) = Uy tan’ + x (26) | 
The variational quantities are related by 
x; = Xe 
Xo = x3+ Kacy) (27) 
ts = —w'x, + Kay 
where K,, is 2.2, o is 25 sec., and wis 0.2 rad./sec. The 


controlled element is also described by the transfer | 


function 


x1(A) 
y(A) 


(28) 


DN ] 
7 K| a + | 
a(A* + ow) 


The results presented demonstrate the following effects: 
(1) variation in initial conditions, (2) mismatch between 
parameters of the controlled member and the estima- 
tion device, and (3) random noise entering the con- 
For all computer solutions, the de- 
The angle I is 


troller directly. 
sired final condition x;p is —50 feet. 
—0.05 rad. The parameters of the estimation device 
are the same and are based on the model of Eq. (27). 
The controlled member and the remainder of the con- 
trol system were synthesized on the D.A.C.L. Gener- 
alized Computer. The integrations in the controlled 
member and the time-varying gains were instrumented 
using rate and position servomechanisms. The con- 
troller was instrumented with a special relay-operated 
device. Hence, the tests used mechanical as well as 
electrical operations in the computation. 

Figure 6(a) demonstrates the effects of variations in 
The length of the control interval 
for each solution is the same, 7 = 20 sec., and the 
desired final condition xp is +7.2 ft./sec. The alti- 
tude condition is met in all cases. The velocity con- 
dition is met for every case except where all three 
Differences in the 
the 


initial conditions. 


initial conditions are positive. 
final conditions would have a 
system performance; although the shape of the flight 
path varies, the final conditions are met unless large 


Then, the performance of the sys- 


similar effect on 


variations occur. 
tem deteriorates when the dynamic limits of the con- 


trolled member are reached. The values of initial 
conditions for Fig. 6(a) are listed in the table of Fig. 
6(b). 


Effects of mismatch between the parameters of the 
controlled member and those of the estimation device 
are illustrated in Fig. 7(a). The variations of the 
parameters used for these solutions are similar to the 
variations that would occur when a particular aircraft 
is operated under a variety of load and climatic con- 
ditions and a variety of mean wind speeds. Because 
variations in the mean wind speed affect the speed of 
the aircraft with respect to the ground, the length of 
the control interval also varies. The parameters of 
the controlled element and the length of the control 
interval are listed in the table of Fig. 7(b). For certain 
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ranges of parameter variations, there is no visible effect 
on the final conditions although the shape of the flight 
path varies. For other ranges of variation studied, the 
control system cannot compensate for the errors in 
estimation caused by mismatch. However, by suit- 
able choice of design parameters, the full range of vari- 
ations that can be tolerated by the control system may 
encompass the variations that are likely to occur in an 
actual controlled member. 

The behavior of the control system when parameter 
mismatch is present is similar to the behavior in other 
practical situations that affect the actuating error in 
such a way that the error can be regarded as the sum 
of two components. The first component is_ the 
actuating error for the terminal control system under 
ideal conditions, and the second is another time-vary- 
ing component. Other situations where the actuating 
error would take this form occur when the prediction 
of the desired response used as a command input varies 
with time, or when approximations are used for the 
time-varying gains in the estimation device. 

Effects of random disturbances are illustrated in Fig. 
8. The length of the control interval for each run is 
20 sec., and the desired final condition x2p is 7.2 ft./sec. 
The random disturbance has a Gaussian amplitude dis- 
tribution with a root-mean-square value of 6 ft. The 
disturbance is obtained by passing nearly white noise 
through a simple-lag filter that has a break frequency of 
100 rad./sec. The dead zone of the controller is 
+().25 ft. Although the root-mean-square value of 
the disturbance that is added directly to the actuating 
error is several times the limit of the dead zone, the 
scatter of final response rates and responses is small. 





(a) Flight paths; (b) 











Fic. 7. Effects of parameter variations. 

Pertinent values: 

A B Cc D* E 

T, sec. 5.5 9 19.4 10.9 8.3 
vep, ft./sec. 16.2 9.2 3.2 7.2 10.2 
Ka 1.6 3.8 1.5 2.2 0.36 
a, sec. 54.5 15.2 20.1 25 92.7 
w, rad./sec. 0.13 0.20 0.25 0.2 0.15 
¢ + 0.06 + 0.08 — 0.06 0 + 0.02 





* No mismatch. 





TERMINAL 


CONTROL SYSTEMS 709 





HORIZONTA TANCE 


Fic. 8. Effects of random disturbances on five flight paths 


Conclusions 


A basic statement of the principles of multicondition 
terminal control has been presented together with a 
possible synthesis procedure for the control equip- 
ment. The method of control is more general than the 
paper explicitly indicates because the principles can be 
extended to include the simultaneous control of several 
linear combinations of response components at a termi- 
nal instant. Control of the particular combination 
used to form the actuating error takes precedence over 
control of the remaining combinations, which are used 
to determine the assumed forcing program. 

The method of multicondition terminal control is 
applicable to any constant-coefficient or time-varying 
linear controlled member. In addition, members with 
input limitations imposed by saturation or undesirable 
response effects, such as aerodynamic stall, can be con- 
trolled. 

The experimental results indicate that moderate 
amounts of parameter mismatch, of variations in initial 
or final conditions, and of random disturbance do not 
deteriorate the terminal performance. 


Appendix I—General Solution to a Set of Linear 
Differential Equations 


The set of linear first-order differential equations 
used as a mathematical model, Eq. (3), is placed in 
matrix form, and the general solution to such a set of 
equations is obtained. The concepts of addition, multi- 
plication, inversion, differentiation, and integration of 
a matrix are used. A description of these concepts can 
be found in the literature.® 
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Fic. 9. The derivative of the magnitude of a quantity. 


The response vector X has as its components the re- 
sponse variables v;, and the matrix 4 is made up of the 
elements a;;._ For a single-input system, the forcing 
vector Y has y as its first component, and the 5; are the 
elements of the first column of B. The matrices -1 and 
B are n X n square matrices and each of the vectors 
X and Y has » components. There are » linearly inde- 
pendent nonzero solution vectors to the homogeneous 
portion of the simultaneous equations’ written in matrix 


form as Eq. (29). 
X AX + BY (29) 
The solution vectors are used to form the columns of a 


matrix ® that also solves the homogeneous portion of 
Eq. (29). 


& = AP (30) 
Postmultiplication of Eq. (380) by ®~' provides an ex- 
plicit expression for 1. 

oh! = 4 (31) 


Differentiation of the identity matrix of order 7 formed 
from © and ®~! yields another relationship for 1, 


—hib-'/dt = A (32) 
because 
1m-! 


d r ; 
@p—!| = hgh! P = = {0 (33) 
dt | | - dt , 0] 


Substitution of the result of Eq. (32) into Eq. (29), 
premultiplication by @~', and rearrangement of terms 
yield 


de 


@-'X + 7 X = &"' BY (34) 
¢ 


The left-hand side of Eq. (34) is a perfect differential, 
and integration of both sides between 7; and 72 and pre- 
multiplication by ®(72) yield 


X(72) = B(r2) B-'(7,) X(7,) + 
P(r.) ®-'(r) B(r) V(r) dr (35 


The component form of Eq. (35) is written as Eq. (5) 
of the paper and is interpreted there. 


Appendix II—Requirements for the Forcing 
Input 


The requirements for the forcing input to be used for 
a nonzero actuating error are derived in this appendix. 
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Principle 3 states: “If the actuating error is not zero 


the control system operates on this error to generate 


a forcing input for the controlled member that will re 
duce the magnitude of the error.’’ For the magnitude 
of a nonzero actuating error x, to decrease, the rate of 
change of xz must be opposite in sign to x, itself, as 
shown by Figs. 9(a) and 9(b). 

The actuating error is defined by Eq. (12). Differ- 
entiation of this expression with respect to ¢ yields 

d d d 


XE = Yipp — Via 


at” dt dt 


(2 
(30 


For the magnitude of x, to decrease, the condition 


d d ae 
sgn Xipp — X14] = —SQN Xz (o4 
dt dt 


must be met. The sign of the derivative of x, is af- 
fected by the forcing input only through the action of 
this input on the rate of change of the estimated re- 
sponse x;4. The rate of change of the predicted input 
X1pp is determined by the input data. 

If the forcing input is so chosen that 

sen — Xi4 = SEN Xz, Xz ~ 0 (38) 

the control system operates to reduce the magnitude 
of the actuating error x, whenever this is consistent 
with amplitude limitations that may be placed on the 
forcing input y. The truth of this statement is demon- 
strated in the next paragraphs. 

The sign of a difference between two quantities is 
determined by the sign of the quantity with the larger 
magnitude. Hence, because 


d d d d d 
a ee aoe” lee a 


the sign of the derivative of the actuating error is con- 
trolled by the sign of the derivative of x, if 


d > d (40 
p Xipp ‘ 
dt = dt _ 
In thi ‘ = v1 (41 
n this case, sgn Xe sgn “VIA 1) 
dt dt 


But, because the sign of (d/dt)x,, is chosen, according 
to Eq. (38), to be the same as the sign of xz, the condi- 
tion of Eq. (37) is met, and the magnitude of x decreases. 

Equation (39) also shows that the sign of (d/dt)x, 
is determined by the sign of the derivative of the pre- 
dicted input, (d/dt)x1 pp, if 


d m d (49 
Xipp X14 a 
dt dt 
For the case where 
d , 
sgn — Xipp = —Sgn Xz (43) 
dt 


the condition of Eq. (37) is still met, and the magni- 











Di 


tit 
lin 
If 

th: 


gerc 
sch: 


9 


Mu 
Jan 


24 
flect 


fron 





1) 


if - 


of 


le 
it 
le 


iS 











MULTICONDITION TERMINAL CONTROL SYSTEMS 711 


tude of xg decreases. When the sign of (d/dt)x;pp is 
the same as the sign of xz, the condition is not met and 
the magnitude of x, increases, contrary to the desired 
action of the control system. This undesired effect is 
lessened, however, by the proper choice of the sign of 
(d/dt)x:4. For the choice of Eq. (38), the magnitude 
of (d/dt) xg is minimized because both (d/dt) x; pp and 
(d/dt) x14 have the same sign. Also, for the desirable 
behavior occurring when Eqs. (42) and (43) hold, the 
magnitude of (d/dt) xg is maximized for the same reason. 

The preceding discussion illustrates the advantages 
of maximizing the magnitude of (d/dt) x4, the rate of 
change of the estimated response. This purpose is 
accomplished by a proper choice of the forcing input. 
Utilizing the results of Appendix I, the estimated re- 
sponse X;4(7\f) can be expressed as 


t 
Via(7,t) = x14(7,0) + f hy(7,r) y(r) dr + 
0 
7 
hyi(T,r) ya(r) dr (44) 
Jt 
Differentiation with respect to ¢ yields 


: Ma(T) = b(T OW — val(d)] (45) 
dt y 
The sign of (d/dt)x,4 is controlled by y(t) if this quan- 
tity is larger in magnitude than y,(¢). If no amplitude 
limit is present, this condition always can be enforced. 
If the amplitude is limited, a forcing program is used 
that has amplitudes less than the limit level. If 


Sublayer Theory (Continued from page 678) 


y(t)| = Vz Xe! #0 (46) 
the sign of (d/dt) x4 is still controlled by yi‘). The 
magnitude of (d/dt) x14 is maximized if the differ- 
ence |y(t) — y,4(t)| is maximized. These conclusions 
lead to the requirements of Eq. (13) for the forcing y(t), 


y(t) | =] lyalt) 


sgn y(t) = [sgn /,(7)t)| [sgn xp(t)] 
y(t)! — (y4(t) | is maximized 
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Relation Between Transition Location 
and Pressure in Separated Shear Flows* 
J. Leith Potter and Jack D. Whitfield 

von Karman Gas Dynamics Facility, ARO, Inc., 
Tullahoma, Tenn. 

April 11, 1960 


T= PRESENT AUTHORS have shown that base-pressure ratio 

(pPo/P.) is essentially a function of transition location when 
transition occurs in the separated wake of similar bodies at equal 
Mach numbers.!:? 
afterbody shapes are similar and local pressures and Mach num 
The hypothetical 


(The bodies need not be entirely similar if 


bers just upstream of the bases are equal.) 
length XY, (Fig. 1) was used previously to correlate base-pressure 
data because of difficulty in determining the mean station of 
transition in a separated shear layer from schlieren photographs, 
particularly under high-Mach-number, low-density conditions 

An example of the correlation of base-pressure data corre 
sponding to transition in the wake between the points of separa 
tion at the base and wake closure downstream is displayed in 
Similar correlations have established for Mach 
The correlation clearly demonstrates that 
However, 


Fig. 2. been 
numbers of 2 and 3.2 
transition location is the dominant factor in this case 
one would like to know how Xy, and X;; are related, since X¢,; is 
the direct variable. 

Inspection of numerous schlieren photographs of the wake 
flows from the ogive-cylinder models used to establish the correla 
tion has revealed the approximate relation between X and X;; at 
free-stream Mach numbers of 3 and 4. As expected, replacing 
X»/D with X,;/D in Fig. 2 simply compresses the distance scale 
by a numerical factor, thereby supporting the correlation pro 


cedure of references | and 2. 
* Sponsored by Arnold Engineering Development Center, Air Research 


and Development Command, USAF, under Contract No. AF 40(600)-700 
S/A 13(59-1) 





Contributors Please Note 


Readers’ Forum items must be confined to the 
equivalent of one page in the Journal. Contribu- 
tions that exceed this limit will be returned to the 
authors for condensation and rewrite. To avoid incon- 
venience and delay, please adhere to the following 
specifications: 

(1) Five, double-spaced, typewritten, manu- 

script pages (8% by 11 in.), including 
wide margins, formulas, and headings, 
equal one printed page. 
For every illustration, deduct at least 
one-half or as much as one manuscript 
page, depending on the size of the 
illustration. 


(2 


— 

















Th 


ber 


figure 
It 
comp 
CaSCS 
to se] 
heren 
trans 
prese 
Furtl 
in th 
numt 
Final 
of Re 
const 
Fre 
laver: 
with 
VW 
recov 


temp 


b oo 
oO 


BASE PRESSURE RATIO, p. /p 
oO 


/p 
oo 


b 


BASE PRESSURE RATIO,p 
oO 














The ratio X,;/X», was found to differ at the two Mach num- 
bers and this suggested the preparation of Fig. 3. The latter 
figure confirms a trend mentioned in reference 3 

It is emphasized that Fig. 3 does not compare bounded and 
completely free shear layers—i.e., the boundary layers in both 
cases were bounded over the distance from the stagnation point 
to separation. Another obvious factor in these data is the in- 
herent variation of pressure along the wake. This means that 
transition falls in regions of varying pressure gradient, but the 
present data do not reveal a systematic effect of this variable 
Furthermore, a supporting sting of diameter 0.3D was present 
in the tests of axisymmetric bodies. Also, the lower Mach 
number data on Fig. 3 are from tests of 2-dimensional bodies 
Finally, the ratio Y,,;/Xy is presented on the basis of equal values 
of Re, for each X,; and Xy ratioed. However, Re; is not 
constant at each Mach number in Fig. 3 

From this evidence, it would appear that separated shear 
layers become relatively more stable under conditions associated 
with higher Mach numbers. Cited data indicate Xi/Xy», = 

\f) to a first approximation when walls are close to adiabatic 
recovery temperature. It is clear from reference 4+ that wall 


temperature can produce a strong secondary effect on this ratio 
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The Effect of Thermal Radiation on the I nviscid 
Hypersonic Flow Over a Blunt Body 


G. A. Bird 

Department of Aeronautical Engineering, The University of Sydney 
Sydney, Australia 

April 11, 1960 


Ww CALCULATING THE DETAILS of the flow over blunt 
bodies at supersonic speeds, it is usual to assume that the 
flow behind the shock wave is adiabatic and that the entropy is 
constant along each individual streamline. However, at hyper 
sonic speeds, the gas behind the shock is heated to extreme tem 
peratures and an appreciable amount of heat is lost from it due 
to thermal radiation 

A useful parameter for assessing the importance of the radicnt 
heat loss is the ratio ® of the heat lost by an element of gas near 
the surface between the stagnation point and the sonic point to 
the change in enthalpy of the element, due to the acceleration of 
gas over this interval. In order to evaluate the order of magnitude 
of this parameter, the temperature 7, density p, sound velocity 
a, specific-heat ratio y, gas constant R, and emissivity € per unit 
layer thickness of the gas will be assumed constant at certain 
average values over the interval. The heat loss per unit mass of 
gas due to radiation is given by 


AQ = (eo]4/p)At 
where o is the Stefan-Boltzmann constant and Af is the time 
taken for the element to traverse the distance between the stag- 
nation and sonic points on the body. The velocity distribution 
between these points is not very different from linear, so to a first 
approximation is 
At = 2Il/a = 2I(yRT)~? 
The emissivity € is given by Rose! and, over the range of interest 
T = 6 — 8,000°K., p = 1 atm.), it is related to some reference 
value by the approximate empirical relationship 
€ = éret( T/T ret) *( p/ Pret 
The change in enthalpy due to the acceleration of the gas is 
Ah = a?/2 = yRT/2 
Therefore, ® = (AQ/Ah) = [(4eere¢ 1T%5)/(yR)*!* prot Tret*] 


and, for air as a real gas at the conditions of interest behind the 
shock, the ratio of the temperature 7 to the ambient temperature 
T, is very nearly a linear function of Mach number 1. 


Therefore, ® = const. (]M*57,6-) 


where] may be regarded as a typical dimension of the body. 
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near the stagnation streamline. 


The author’s? numerical method of solution of blunt-body 
flow fields was modified to allow for this heat loss by the addition 
of the appropriate term to the one-dimensional energy equation 
which is applied along each stream tube and the calculation of the 
progressive change in entropy along each stream tube. The 
method was then used to calculate the details of the flow behind 
a paraboloidal shock wave of nose radius 5 ft. traveling at a 
Mach number of 20 in air at 220°K. and 0.01 atm. Fig. 1 shows 
the path on a thermodynamic chart of an element of gas very 
close to the stagnation streamline (minimum velocity 25 ft./sec.) 
when thermal radiation is taken into account and also when it is 
neglected; and it is seen that the parameter ® is here equal to 
2/3. The greater part of the heat loss and consequent change in 
entropy occurs in the stagnation region where the flow velocity is 
small. However, there was no detectable change in the geometry 
of the flow or in the pressure distribution on the body, and the 
only significant change was in the temperature of the gas layers 
near the surface. 

The near constancy of the flow geometry and pressure distribu- 
tion can be explained by the following argument. The flow is 
constructed as a series of one-dimensional stream tubes and the 
geometry of the flow remains unchanged unless there is a change 
in the flow velocity or density along them. The radiant heat loss 
causes a Significant drop in temperature only in the stream tubes 
which pass near the stagnation point (such as that in Fig. 1). 
The small change in area of these stream tubes due to the re- 
sultant density change does not produce any significant change in 
the overall flow pattern. Therefore, as long as there is no change 
in the flow velocity, the pressure gradients and pressures will re- 
main almost unchanged; and it can be easily shown® that, ir- 
respective of the amount of heat lost, there is no change in veloc- 
ity unless there is a change in pressure. 

This study indicates that, for practical hypersonic blunt-body 
flows, the velocity gradient outside the boundary layer will not 
be significantly altered by thermal radiation effects, and the 
effects on the temperature can be calculated by a step-by-step 
construction along the stagnation streamline using the pressure 
distribution given by theories which neglect radiation. These 
conclusions will only be affected if the parameter ® is increased 
by at least an order of magnitude above its value in the present 


example. 
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A Note for Finding Complex Roots 
of a Special Kind of Quartic Equation 


F. A. Lee 
Structural Dynamics Group, Structural Research Unit, 
Boeing Airplane Company, Transport Division, Renton, Wash 


April 11, 1960 


’ j ‘WE ALGEBRAIC EQUATION of the following type 


xt + Aox? + Ayx + Ay = O (1) 


has been found to arise in the flutter analysis of low-aspect-ratio 
wings! and supersonic panel flutter of a evlindrical shell of finite 
length.2 In the references, the equation was solved by a trial 
and-error method. Here we propose to present an explicit 
method for finding the complex roots of this kind of quartic equa 
tion. 

Assume Eq. (1) can be factored into 


(x? + dx + by )(x? + aux + ao) = (2 


By equating the coefficients of different powers of x of Eqs. (1 


and (2), we have 


a, +b = 0 (3) 
a,b; tT do T by = wl, (4) 
ayby + aby = Ay (5) 
doby = Av (6) 
From Eq. (3), let a; = —),; = a, then Eqs. (4) and (5) can be re 
duced to 
ay + by = A» T a‘ (4 
— iy + hy = A, a (S 
Therefore, 
2a) = A» +a*— (A, a (9) 
2b) = Ao + a? + (A, a} (10) 


Combining with Eq. (4), we get the governing equation for a? as 
follows: 

4Ay = (Az + a?)? — (A;7/a? 
This equation can be written alternatively as 


B® + 2A, B? + (Az? — 4A40)B — Ai? = 0 


where 


The solution of the above cubic equation is well known and the 


solutions are 


Bb = A+B — (p/3) / 

Bo = —|(A + B)/2] + [(A B)/2|~ -3 p/3) > (il 

8; = —((A + B)/2] — ((A — B)/2] V—3 — (p/3)) 
where 


3 

A = V —(b/2) + V(b2/4) + (43/27 
3/ 

B = V —(6/2) + V (b2/4) + (03/27 


(12 


a = (1/3)(38q — p?) 
b = (1/27)(2p3 — 9pq + 27r) 


It is noted here that 


(i) if (b?/4) + (a%/27)> 0 there will be one real root and 
two complex-conjugate roots. 

(ii) if (b2/4) + (a3/27) = 0 there will be three real roots of 
which two at least are equal. 

(iii) if (b?/4) + (a3/27) < 0 there will be three real and 
unequal roots. 


Only real and positive roots are of interest here. 


Once the values of 8 are determined, the values of ap and 4p can 
be found by Eas. (9) and (10). Therefore, we have the original 
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equations exactly separated into two quadratic equations, for 
which the roots are readily available by means of the quadratic 
formula 

rhe following example is a simple illustration of the above 


method. Let us consider the equation 


‘io+t+xt+x—1=0 
From the relations given under I 12 
i,=1, A, =!1 1, = l 
p =2, q=5, r= 1 
a = (1/3)(3 & 5 — 4) = 3.6667 
b =(1/27)[2x* 8 -9xX2x*5+4+27(-1)] = —3.7407 


Applying criterion (i) at top of page 
(62/4) + (a3/27) = 5.3232> 0 
Therefore, we find that there will be one real root and two com- 


plex-conjugate roots To calculate the roots 


3 
A = V 1.8704 + 2.3065 = 1.6105 
3 
B = V 1.8704 — 2.3065 = —0.7583 
~-pB=0 1855, a = 0.4307 


By Egs. (9) and (10) 


a = —0.5682 
bo = 1.7537 
Eq. (2) can be written 
(x? + 0.4807 x — 0.5682)(x? — 0.43807 x + 1.75387) = 0 


The roots of the above quadratics are 


x, = 0.5689 
vm = —0.9996 
v; = 0.2154 + 1.30167 


vy; = 0.2154 — 1.30161 
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2 Holt, Maurice, and Strack, Sainsbury L., Supersonic Panel Flutter of A 
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A Particular Solution to the Turbulent- 
Boundary-Layer Equations 


George K. Walker 

Aerophysics Engineering Operation, MSV Department, 
General Electric Company, Philadelphia 4, Pa 

April 18, 1960 


INTRODUCTION 
ECENTLY, several investigators have extended the turbulent, 
compressible, zero-pressure-gradient boundary-layer rela- 
tions to include the case of a blunt body with a finite pressure 
gradient. Most investigators employ the momentum and/or 
energy integral equations, the semi-empirical Blasius flat-plate 
skin-friction coefficient with a modification for compressibility 


effects, and several of the following assumptions: 


(a) Unit Prandtl number 

(b) Incompressibility (p/pe = 1) 

(c) A power-law velocity profile 

(d) A Crocco-type enthalpy-velocity relation 
(e) A highly cooled wall 

(f) Reynolds’ analogy 


The present derivation replaces assumptions (a) to (f) by a single 
relation which greatly simplifies the solution of the equations 


ANALYSIS 
The energy integral equation is given! as follows: 
- 


! 
C/es(o UhOr oe) = rq (1) 


FORUM 715 


where ¢ is the enthalpy thickness, defined as 


a - 
Cs pil ( 1 - h \ dy . 
0 é , 


where the superscript 0 denotes totalenthalpy. Also k 0 fora 


two-dimensional body and k 1 for a body of revolution 
The relationship between heat transfer and skin friction is 


assumed to be described by the function G, where 
G(s )(e¢/2)(h, — h 3) 


q = (pe U, Pr 


t transfer and 


where G(s) is defined as the departure of the he: 


skin friction relation from its zero-pressure-gradient value 
(G = i) 

Now the Eckert zero-pressure-gradient form of the skin-friction 
coefficient is assumed adequate for the case of finite pressure 
gradient; therefore, 

Cr/2 = 0.013(pe/ pe U0): e!- 5 = N(s)/O°-? (4 


where 6 is the usual momentum thickness, and where 


S v-5 (5 


€ (a” Me )-2( o*/ pe} 
The starred quantities indicate properties that are to be evaluated 
at the local values of pressure and reference-enthalpy 
The relationship between ¢ and @ is assumed to be described by 
the function ¢, where 


¢ = [Chr — hy)/(Prv*he) \e(s)e (6) 


where ¢(s) is defined as the departure of the ¢ and @ relation from 
its zero-pressure-gradient value (¢ = 1) 
When Eggs. (3), (4), and (6) are substituted into Eq. (1), in 


tegration vields 


ec. { s 
[1 25 a [peUer*(hy — hw) )®NY-%Gds 


ray 
pel er(h, — hy) 


A convenient solution to Eq. (7) is obtained for the particular 
solution where 


G = 1/%-% (8) 


The function ¢ can ther be evaluated by solving the momentum 
integral equation for 6, where the skin-friction coefficient is given 
by Eq. (4), and subsequently dividing Eq. (7) by @: 


& oR 
: 7 440-25 ¢1-25( fp, — fp.,,)1-25y1-2kds | 
edi lf pe Ueme®-*e y hee i 
_ U, 0 (9) 
— hy) ’ + 1.95 +2.25,. 0.25.1.25y1. 25%) 
(hy ly i pel p1eRbM +2. 254 0.2561 25y 1. Wied s 


where #/ is the form factor 

Now Eggs. (1), (6), (8) and the momentum-integral equation 
represent four independent equations and five dependent vari- 
ables: gy, 0,G,¢, and H. Therefore, unless additional assump- 
tions are made, H must be given in order to complete the solution 
However, it is interesting to note that substitution of Eqs. (4), 
(7), and (8) into Eq. (3) yields the result that g¢ is independent of 


Ve fe 


1.25/ 1 


0.0296 = peUe[mer™( hr — hw) |? (h. — hy) 


a 
t >»* s 0.2 
Pr [ f. pe U epee (hh, — Ny )'s 257-254 as| 
0 


Therefore, an equation for turbulent heat transfer to a blunt 
body has been derived without resorting to any of the usual 
assumptions, (a) to (f). However, the solution does depend on 


the assumption given in Eq. (8) 


DISCUSSION 
Since the present analysis is similar in some respects to that of 
reference 2, it is of interest to compare Eq. (8) to the correspond- 
ing relation in reference 2, namely 
Ge=14+ [(n+3)/(n +1)] — (1 +3)/(m + Dt (11) 
where » is the exponent of the assumed power-law velocity profile, 


Thus, it is seen that Eqs. (8) and (11) are identical for 
1, and both decrease for 


(n =7) 


the case of zero-pressure gradient, ¢ = 
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increasing ¢. However, the two equations differ in magnitude 
for £ sé 1. 

It should be mentioned that although the solution indicated by 
Eq. (8) is admittedly one of convenience, the solution indicated by 
Eq. (11) is not necessarily more valid. The latter solution is 
largely dependent on the derivative of an assumed polynomial for 
the relationship between total enthalpy and velocity in the 
boundary layer. The use of this mathematical step does not 
always guarantee close agreement between the approximate 
solution and the actual case. Actually, experimental data are 
needed to determine if either Eq. (8) or (11) is valid. 

A sample calculation for a hemisphere has been made for 
typical re-entry conditions, and the results are shown in Figs. 1 


and 2. 


REFERENCES 
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Some Real-Gas Flow Parameters for Air 


Wayne D. Erickson 
Head, Thermal Projects Section, Langley Research Center, NASA, 
Langley Field, Va. 


April 28, 1960 


"T OPERATION of air nozzles at high Mach numbers requires 

that the stagnation air be heated to high temperatures to 
avoid condensation of its constituents upon expansion. In 
some cases it is also desirable to operate at high stagnation 


pressures. 


SCIENCES 1960 


SEPTEMBER, 


Flow-parameter ratios for air as a real gas in chemical equilib- 
rium have been calculated for Mach numbers of 10 and greater 
for stagnation temperatures and pressures up to 4950°R. and 
1000 atm. by using the thermodynamic properties of air pre 
sented in reference 1 and extensions thereof. (The data of refer 
ence 2 support the validity of the extended properties.) This 
note is concerned with but three flow ratios, viz., (1) free-stream 
static temperature to stagnation temperature, (2) free-stream 
static pressure to stagnation pressure, and (3) total pressure be 
hind a normal shock to stagnation pressure, defined, respectively, 
as, /’2,and 73. For areal gas, each of these ratios is in general a 
function of both stagnation temperature and pressure as well as 
Mach number. This three-fold dependency presents the awk 
ward problem of presenting the flow ratios for a real gas once they 
are calculated. 


Reference 3 presents, in addition to the tabulated ideal-gas 
flow ratios for constant heat capacity, which are of course a single 
function of Mach number, correction factors for the various flow 
ratios based on calculations which include only the variation of 
heat capacity with temperature. These correction factors—i.e., 
the quotient of the ideal-gas flow ratio as calculated with regard 
for variable heat capacity and the ideal-gas flow ratio as calcu 
lated with constant heat capacity—are plotted for several stag 
These 


f the 


nation temperatures as a function of Mach number 


plots show that the correction factors are independent « 
Mach number for Mach numbers greater than about 7 


With the anticipation of this independency of Mach number 
for a real gas, the calculated real-gas flow ratios for various 
Mach numbers greater than 10 for a given stagnation tempera 
ture and pressure were divided by the corresponding ideal-gas 
ratios tabulated in reference 3. The resulting quotients (correc 
tion factors, defined as g = r/r’, where r’ is a particular flow ratio 
calculated for an ideal gas with constant heat capacity) for the 
three above-mentioned flow ratios were found to be independent 
of Mach number and could, therefore, be presented for various 
stagnation temperatures as a function of stagnation pressure on 
a single plot. These plots of correction factor against stagnation 
pressure for various stagnation temperatures shows that the 
correction factor is very nearly a linear function of stagnation 
pressure for a given stagnation temperature. For a particular 


flow ratio, the correction factor can, therefore, be represented as 
¢ = go tmp, (1) 


where ¢, and m are functions of stagnation temperature alone and 
Pt. is the stagnation pressure in atm 
Now when the values of g for a given ratio were plotted 
against the stagnation temperature, it was found to be linear in 
the stagnation temperature and could therefore be represented as 
nA 2 > ult (2) 
where a and uw are constants for a given flow ratio and 7;,; is the 
stagnation temperature in °R. The constants are tabulated in 
Table 1. 


TABLE | 
Correction Factor a uw X 104 
¢) 0.948 0.516 
¢2 1.181 —1.14 
¢g 1.188 —1.13 
TABLE 2 
Correction Factor BX 10! n X 108 
¢1, (T1.1 < 4590°R.) 0.958 —1.185 
2, (T1.1 < 4080°R.) 4.196 —7.67 
(Ti. > 4080°R.) 2.201 —2.76 
¢3, (T1.1 < 4230°R.) 3.887 —6.76 
(T1,; > 4230°R.) 2.563 —3.66 
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A plot of m against stagnation temperature for a given ratio did 
not yield a relationship that could be expressed analytically but 
by segmenting the plot into two regions, a linear approximation 
was found in each region between m and the stagnation tempera- 
ture. For each temperature range then 


m=B+nT1, (; 


w 


where 6 and y are constants and shown in Table 2. 

The region in which these equations apply is shown by the 
shaded area in Fig. 1. Also shown in Fig. 1 by the nearly vertical 
lines are approximate combinations of stagnation pressure and 
temperature at which condensation for equilibrium conditions of 
dry air will just begin for various free-stream Mach numbers. 
To the right of these lines, condensation-free flow will exist for a 
respective Mach number. 

The procedure for finding the three flow-parameter ratios con- 
sidered herein is best shown by an example. Consider a case 
where 7;,, = 3510°R., pr, = 400 atm. and the measured pitot 
total pressure is 0.145 atm. The ratio 7; is 0.145/400 or 3.62 
10-4. The correct values from Table 1 and Table 2 for deter- 
mining ¢; are a = 1.188, 4 = —1.13 XK 1074,8 3.887 =X 10-+, 
and » = —6.76 X 1078. Therefore, from Eq. (2) 


gos = 1.188 — 1.13 K 1074 X 3510 = 0.791 
and from Eq. (3) 
m3; = 3.887 X 107-4 — 6.76 K 107-8 K 3510 = 1.515 K 10~4 
From Eq. (1) 
g: = 0.791 + 1.515 KX 10-4 K 400 = 0.852 
The ideal-gas value of r3’ is then 
3.62 X 10~—4/0.852 or 4.25 X 10-4 


The corresponding value from reference 3 for J/, is therefore 
15.10, which is the Mach number based on the real-gas calcula- 
tion. The other two flow ratios may be found likewise. 

The flow parameters calculated by this scheme were found to 
differ no more than +0.5 per cent from the values determined 
directly from the thermodynamic data for air as for the limited 
stagnation conditions shown in Fig. 1. However, in most cases, 
the difference is less than +0.2 per cent. 


REFERENCES 

1 Hilsenrath, Joseph, and Beckett, Charles W., ¢ al., Tables of Thermal 
Properties of Gases, NBS Cir. 564, U. S. Dept. Commerce, 1955 

2 Hall, N. A., and Ibele, W. E., Thermodynamic Properties of Air, Nitrogen 
and Oxygen as Imperfect Gases, Tech. Paper No. 85, Eng. Exp. Station, Univ 
Minn., December, 1951. 

3 Ames Research Staff, Equations, Tables, and Charts for Compressible 
Flow, NACA Rept. 1135, 1953 


FORUM 717 


Axisymmetric Bending Stresses in Solid 
Circular Plates With Thermal Gradients* 


Marvin Forray and Malcolm Newman 


Design Specialist and Principal Engineer, Structures 
Applied Research and Development, Republic 
Aviation Corporation, Farmingdale, L.I., N.Y. 


May 20, 1960 


bon PROBLEM of the axisymmetrical flexural stresses in 
clamped and simply supported solid circular plates due to 
thermal gradients was considered in a previous paper.!' How 
ever, the frontal attack of solving the partial differential equa 
tion for the deflection was not employed. It is the purpose of 
this note to solve the problem directly, and to present curves and 
formulas which will be useful to the practicing engineer. While 
the discussion will be restricted to the axisymmetric problem, it 
should be noted that the asymmetrical case may be treated by 
the procedure in references 2 and 3 

A diametrical section of the plate is shown in Fig. 1. The 
equilibrium deflection equation for an elastic plate subject to an 
axisymmetric temperature distribution, which varies linearly 
through the thickness, takes the form (reference 4 


V tw {( 1 + v)a/h IV?7 p (1) 
where 7p is the temperature difference between the upper and 


lower faces, respectively, v is Poisson’s ratio, a is the linear coeffi 


cient of expansion, and in polar coordinates 


ws (< ‘ ] “\(< ld 
~~ \dr? r dr} \dr? ? =) 


It is assumed that the temperature is expressible by a con 


vergent power-series expansion of the form 


* This research was supported by the USAF under Contract No. AF33 
616)-6066 monitored by Mr. C. Richard, WCLSS-T30, Wright Air Develop 


ment Division 
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ox 


axrk +C 


T(r,s) = 
h kK=0 
(ax, C are constants) so that the temperature difference Tp is 
given by 
Tp = = axr® (2) 
K=0 
The solution of Eqs. (1) and (2) is 
Cir? ; , (1+ v)a = axr® *2 
w= rm + C. + = : (3) 
a* h K=0 (AK + 2)? 


where C; and C, are arbitrary constants which must be deter- 


mined by the boundary conditions. 
(A) CLAMPED PLATE 


The solution of Eq. (3) for the clamped plate [w = dw/dr = 0 


(1 + via aK rr 5 aa a 2 
w= = ~ (rE ** + KaX** — (K + 2)r*e*] 
Ph K-0 (K + 2)? 
(4) 


The bending moments in the radial and transverse directions 
M,, Mg are given by (reference 4) 


— Eh?a 
M, = 
121 — » 
x ; = : K 7 
 B _ Jc +140“) ~a+nb—t (5) 
rok +24 a f a 
— Eh*a 
Mo 
1241 — » 
axa® ) s (: x \ _ . 
4 [1 + o(K + 1) —(l+v)7—Tbd (6) 
K-0K +2) . MI ) " "9 





The deflections and moments in nondimensional form when 


Tp is expressed by the monomial Tp = axr* are 


wh (1 + v) r\ +3 K (: ? K 
—- = — -_ + 1 : = aie (7) 
aaT p(a) (K + 2)? a 2 a 2 
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NN, 1 (: K 1+» 
+ (S 
Eh?aT p(a)/12 K+2 a l— p 
M l r\K 1 + 
ae = — cen (! 4. “| «9 
Eh? aT p(a)/12 K+2 a l— vp 


where 7/)(a) is the temperature difference at r = a 
nondimensional deflections and moments are presented in Figs. | 


and Ztor K = 0,1,2,...5:;7 = 03 
be used for 7p given by polynomials in r 


Superposition may then 


(B) SIMPLY SUPPORTED OR FREE PLATE 
For the simply supported or free plate subject to a monomial 
temperature distribution, add to the solution for clamped plate 
Eqs. (7)-(9) the following: 


wh (r/a)?> — 1 . 
rye: = = (relative to boundary ) (10) 
a’aT p(a) K+2 
M, Mo 2 
= (11 


Eh®aT p(a)/12 Eh? aT p(a)/12 (K + 2)(v — 1) 


Example 


Given: Aclamped aluminum plate, with 


a= 10in., h = 0.20in., E = 107 psi, 
a = 12 X 10-%in./in./°F., v = 0.3 


where the radial variation of temperature difference Tp in °F. 
is given by 7p = 200 + (r?/2). 

Required: Maximum deflection and stress. 

Solution: If the plate is clamped and the coefficients of the 
polynomial expressing the temperature difference are all of the 
same sign, as above, then the maximum magnitude of deflection 
occurs at the center of the plate while maximum magnitudes of 
radial and tangential moments occur at the boundary. Let 
and ” denote the contribution of the two terms in Tp 

From Figs. 1 and 2(A = 0, K = 2, respectively), 

, M,’ Mo’ 
w' =0, oer = ae = 1.43, 
(Eh? aT p'(a)/12) [Eh?aT p'(a)/12] 
w"h(K + 2)? M,” " 
= 1.30, ge = 0.71, 
rja=0 [Eh?aT p"(a) 12] r/a=1 


ara T pn" ( a) 


Mo” 1.91 
[Eh?aT p"(a ) /12] rja=1 sid 


Substituting Tp’(a) = 200, Tp"(a) = 50, and the given data 
yield 
M,'|r/a=1 Mo'\ry,a=1 = 114.4 in.-lb. /in. 


w”|rJa=o = 0.024 in., M,” = 14.2 in.-lb./in., 


so that finally 
Wirja=o = [w’ + w"Jr/ea=0 = 0.024 in., 


Mr\yja=1 = [M,’ + Me" )y/a=1 = 128.6 in.-lb./in. 


Mo)\r/a=1 = [Mo’ + Mo" ]+/e=1 = 188.6 in.-Ib./in 
6 zs , 

=— | = 20,790 psi 
h? r/a=1 
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READERSY’ 


Boundary-Layer Transition Measurements 
at Mach Numbers from 5.4 to 7.4 


Fred Landis,* Martin R. Fink,** and Murray H. Rosenberg*** 
United Aircraft Corporation, East Hartford, Conn. 
April 11, 1960 


A ttnovce extensive supersonic measurements of transition 

Reynolds numbers are available from tests of two-dimen- 
sional adiabatic models below a Mach number of 4.5, only one 
set of data at a higher Mach number (.7 = 5.8) has been pub- 
lished (reference 1 These data indicate an increase of transi- 
tion Reynolds number with increased Mach number for adiabatic 
wall conditions. In order to determine whether this trend ap- 
plies at hypersonic Mach numbers, tests on a two-dimensional 
model (Fig. 1) were conducted in the United Aircraft Corporation 
Research Laboratories 6 & 6-in. hypersonic continuous-flow tun- 
nel, under the sponsorship of the Pratt & Whitney Aircraft 
Division. The upper and lower surfaces of the model were in- 
strumented with thermocouples to determine the axial distribu- 
tion of recovery factor and also with static pressure taps to meas- 
ure the pressure distribution on the model. The model surfaces 
were ground and polished to approximately 3-5 microinch rms 
roughness and showed no transverse scratches under visual in- 
spection. The Mach number at the transition location was 
varied from 5.4 to 7.4, and the unit Reynolds numbers at the 
transition location varied from 0.2 & 10°/in. to 0.6 X 10°/in. 
Ali tests were conducted at essentially adiabatic wall conditions 
by allowing sufficient running time to reach thermal equilib- 
rium. Transverse contamination from the turbulent sidewall 
boundary layers did not reach the instrumented region of the 
model. The minimum transition distance was defined as the dis- 
tance from the leading edge to the point where the temperature 
recovery factor began to increase with axial distance. 

The flow properties adjacent to the boundary layer varied along 
the model because of a Mach number gradient in the wind tunnel 
and because of curvature of the forward portion of the model. 
Calculations showed that the expansion waves generated by the 
tunnel contour were essentially unaffected by the weak bow wave 
generated by the model. The flow turning angles produced by 
these expansion waves were added to the physical turning angle 
of the model, and inviscid pressure distributions were computed 
by the first-order shock-expansion method. Following the sug- 
gestions of reference 2, the calculated effects of leading-edge 
bluntness and viscous interaction were linearly superimposed 


upon these distributions. Satisfactory agreement was obtained 


* Consultant, Pratt & Whitney Aircraft Division, also Associate Pro 
fessor of Mechanical Engineering, New York University 
** Supervisor, Missile Aerodynamics Group, United Aircraft Corporation 


Research Laboratories 
*** Research Engineer, United Aircraft Corporation Research Labora 
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Fic. 1. Two-dimensional model. 
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Fic. 2. Comparison of calculated and measured static pressure 
distributions 
8.0 -- — ; 
l k 
sts REYNOL! ii — TRANSITION | 
6.0 SYMBOL REF NO. FINCA = aa * 5 s ro eascaress 
‘ o 1 300,000 r } 
° 4 200,000 
5.0}, + 4 —p-a——___ - 
4.0-—_—_+-- | . ts 
3.0-———+- | . aft 4 
| 





UAC DATA 

SOLID UNFLAGGED, a=0* | 
SOLID FLAGGED, a=2.17° 
REYNOLDS 
SYMBOL NO. /INCH 
200,000 
300,000 
409,000 
500,000 


___ 600,000 


<ae@rpace 





TRANSITION REYNOLDS NUMBER x 10° 


5 
MACH NUMBER 


Fic. 3. Transition Reynolds numbers 


between calculated and measured static pressure distributions 
At all Mach 


numbers, the static pressure gradients on the model were small 


(Fig. 2), thus verifying the calculation procedure 
in the axial region in which transition occurred. As indicated in 
reference 3, the effects of pressure gradient on boundary laver 
stability decrease with increasing supersonic Mach number and, 
therefore, the effects of these pressure gradients on transition 
Reynolds number were assumed to be negligible 

The measured transition Reynolds numbers are presented in 
Fig. 3 and compared with the empirical curve given in Fig. 8 of 
reference 4. The solid unflagged symbols represent the results of 
tests at an angle of attack of zero degree and the solid flagged 
symbols indicate data obtained at an angle of attack of 
2.17°. The small difference between both sets of data indicates 
that the slight favorable pressure gradient just upstream of the 
transition region (axial distances of 10-15 in. in Fig. 2) did not 
materially affect the results. Where transition did not occur on 
the model, Reynolds numbers based on model length are shown in 
Fig. 3 with vertical arrows to indicate conservative lower design 
limits. All transitions shown in this figure occurred upstream of 
the region affected by the reflected bow shock wave 

As shown in Fig. 3, the transition Reynolds number increased 
greatly with an increase in Mach number, and increased with an 
increase in unit Reynolds number at constant Mach number 
The transition Reynolds number reported in reference 1 is about 
20 per cent higher than would be expected from extrapolation of 
these data to the same unit Reynolds number. Although this dis- 
crepancy may be due to experimental inaccuracies, it might also 
occur because the model described in reference 1 was adjacent to 
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a partially laminar wind tunnel nozzle boundary layer. The 
nozzle wall boundary layer was turbulent in the tests reported 
herein, and thus the destabilizing effects of tunnel noise (reference 


5) would be more pronounced. 


CONCLUSION 
The following conclusions can be obtained from these test re- 
sults: 


(1) The transition Reynolds number increases with increasing 
Mach number for smooth, two-dimensional, essentially flat 
plate models at Mach numbers above 5. 

(2) In wind tunnel tests, the transition Reynolds number in- 
creases with increasing unit Reynolds number. 

(3) Small favorable pressure gradients upstream of the transition 
location do not appear to alter the transition location signi 


ficantly at hypersonic speeds. 
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Dimensionless Parameters for Viscous 
Similarity* 
Andrew G. Hammitt 


Senior Research Associate, Gas Dynamics Laboratory, 
Princeton University, Princeton, N.J. 


April 29, 1960 


Te REYNOLDS NUMBER is the well known parameter for 

determining the importance of viscous effects in any flow 
problem. It is generally assumed that the Reynolds number, 
calculated for any characteristic length and any typical flow 
condition, is sufficient to define the viscous scale of the problem. 
However, this is not generally the case at hypersonic Mach 


numbers. 

In general, dimensional analysis tells us that for a perfect gas, 
dimensional similarity will be governed by geometry, 7, y, Re 
and Pr. These similarity parameters are sufficient as long as the 
heat-transfer and viscous parameters can be expressed as a single 
dimensional parameter. If more than one dimensional param- 
eter is required to describe the viscous or heat-conduction effects, 
then more dimensionless parameters will be needed to ensure 
similarity. 

* This work was supported by AF Contract 49(638)-465, ARDC, Fluid 
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This point can be illustrated by the use of viscosity alone, but 
the results would be equally valid for heat conduction or any 
other parameter behaving in the same manner. Viscosity of a 
gas is, in general, a function of the temperature and pressure. 
For most practical cases, it may be considered to be a function of 
temperature alone and is often represented as a power law 


= pel T oe ag (1) 


or a Sutherland law 


2 (7 
Me ot er? T)\s 


These are both two-parameter representations of viscosity as a 
function of temperature. The two-parameter equations are less 
limited than a o1e-parameter system, but they are not sufficient 
for all cases. Using a two-parameter viscosity law, dimensional 
analysis shows that two nondimensional parameters are needed 
to describe viscous similarity, Re and 7/.S for a Sutherland-type 
viscosity law. The new dimensionless parameter introduced is 
T/S, which can be important when large temperature ratios 
occur in the flow. 

Consider the example of the hypersonic flow about a blunt body 
and calculate the Reynolds number at (2) behind a normal 
shock as a function of the Reynolds number at (1) ahead of 


Ll + (S/T) (71\" 
Re. = Re ae oe (3) 
1 + (S/7;) \T> 


the shock wave 


If Re, is considered to be one nondimensional viscous similarity 
parameter, S/7, can be taken as the other. If 7,)/72 << 1, then 
S/T, is important. For air, either in a wind tunnel or the 
atmosphere, S/7, & 1 so if 72/7,>>1, S/T2<<1 


l Ty 
Res ial Re, = — 
1 + (S/T) \T? 


Rez depends principally on 1/{1 + (.8/7;)]. 
parameter for air at typical atmospheric and hypersonic wind- 
tunnel static temperatures and for helium at wind-tunnel tem- 


The values of this 


peratures are as follows: 
S(°K.) W4(°K.) S/T, 1/{1 + (S/T 
Air, Atmospheric 190 450 0.42 0.70 
Air, Hypersonic Tunnel 190 65 2s 0.29 
Helium, Hypersonic 
Tunnel 38 5.5 0.15 


Viscous simulation of conditions behind the shock wave may be in 
error by factors of three or more if Re; is held constant. The 
application to helium is questionable because the Sutherland-type 
law is a poor fit to the data over a large temperature ratio and 
similarity between helium and air cannot be obtained since y is 
not the same for the two gases. 

Since Rez is an important parameter in determining the relative 
sizes of the shock-wave thickness, boundary-layer thickness, and 
separation distance, the same Re, may give considerably different 
results depending on the value of S/7) or equivalent parameter. 
In some problems, it may be possible to achieve good simulation 
by matching the particular Re which is most critical, but this 
cannot be a general procedure for obtaining complete similarity. 
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